home *** CD-ROM | disk | FTP | other *** search
/ Enter 2004 January / enter-2004-01.iso / files / maxima-5.9.0.exe / {app} / share / maxima / 5.9.0 / src / ellipt.lisp < prev    next >
Encoding:
Text File  |  2003-02-09  |  102.2 KB  |  3,622 lines
  1. ;;; -*-  Mode: Lisp; Package: Maxima; Syntax: Common-Lisp; Base: 10-*- ;;;;
  2. ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  3. ;;;     The data in this file contains enhancments.                    ;;;;;
  4. ;;;                                                                    ;;;;;
  5. ;;;  Copyright (c) 1984,1987 by William Schelter,University of Texas   ;;;;;
  6. ;;;     All rights reserved                                            ;;;;;
  7. ;;;                                                                    ;;;;;
  8. ;;;  Copyright (c) 2001 by Raymond Toy.  Replaced everything and added ;;;;;
  9. ;;;     support for symbolic manipulation of all 12 Jacobian elliptic  ;;;;;
  10. ;;;     functions and the complete and incomplete elliptic integral    ;;;;;
  11. ;;;     of the first, second and third kinds.                          ;;;;;
  12. ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  13.  
  14. (in-package "MAXIMA")
  15. ;;(macsyma-module ellipt)
  16.  
  17. ;;;
  18. ;;; Jacobian elliptic functions and elliptic integrals.
  19. ;;;
  20. ;;; References:
  21. ;;;
  22. ;;; [1] Abramowitz and Stegun
  23. ;;; [2] Lawden, Elliptic Functions and Applications, Springer-Verlag, 1989
  24. ;;;
  25. ;;; We use the definitions from Abramowitz and Stegun where our
  26. ;;; sn(u,m) is sn(u|m).  That is, the second arg is the parameter,
  27. ;;; instead of the modulus k or modular angle alpha.
  28. ;;;
  29. ;;; Note that m = k^2 and k = sin(alpha).
  30. ;;;
  31.  
  32. ;;
  33. ;; Routines for computing the basic elliptic functions sn, cn, and dn.
  34. ;;
  35. ;;
  36. ;; A&S gives several methods for computing elliptic functions
  37. ;; including the AGM method (16.4) and ascending and descending Landen
  38. ;; transformations (16.12 and 16.14).  We use these latter because
  39. ;; they are actually quite fast, only requiring simple arithmetic and
  40. ;; square roots for the transformation until the last step.  The AGM
  41. ;; requires evaluation of several trignometric functions at each
  42. ;; stage.
  43. ;;
  44. ;; In addition, the Landen transformations are valid for all u and m.
  45. ;; Thus, we can compute the elliptic functions for complex u and
  46. ;; m. (The code below supports this, but we could make it run much
  47. ;; faster if we specialized it to for double-floats.  However, if we
  48. ;; do that, we won't be able to handle the cases where m < 0 or m > 1.
  49. ;; We'll have to handle these specially via A&S 16.10 and 16.11.)
  50. ;;
  51. ;; See A&S 16.12 and 16.14.
  52.  
  53.  
  54.  
  55. (flet ((ascending-transform (u m)
  56.      ;; A&S 16.14.1
  57.      ;;
  58.      ;; Take care in computing this transform.  For the case where
  59.      ;; m is complex, we should compute sqrt(mu1) first as
  60.      ;; (1-sqrt(m))/(1+sqrt(m)), and then square this to get mu1.
  61.      ;; If not, we may choose the wrong branch when computing
  62.      ;; sqrt(mu1).
  63.      (let* ((root-m (lisp:sqrt m))
  64.         (mu (/ (* 4 root-m)
  65.                (lisp:expt (1+ root-m) 2)))
  66.         (root-mu1 (/ (- 1 root-m) (+ 1 root-m)))
  67.         (v (/ u (1+ root-mu1))))
  68.        (values v mu root-mu1)))
  69.        (descending-transform (u m)
  70.      ;; Note: Don't calculate mu first, as given in 16.12.1.  We
  71.      ;; should calculate sqrt(mu) = (1-sqrt(m1)/(1+sqrt(m1)), and
  72.      ;; then compute mu = sqrt(mu)^2.  If we calculate mu first,
  73.      ;; sqrt(mu) loses information when m or m1 is complex.
  74.      (let* ((root-m1 (lisp:sqrt (- 1 m)))
  75.         (root-mu (/ (- 1 root-m1) (+ 1 root-m1)))
  76.         (mu (* root-mu root-mu))
  77.         (v (/ u (1+ root-mu))))
  78.        (values v mu root-mu))))
  79.   (declaim (inline descending-transform ascending-transform))
  80.  
  81.  
  82.   ;; Could use the descending transform, but some of my tests show
  83.   ;; that it has problems with roundoff errors.
  84.   (defun elliptic-dn-ascending (u m)
  85.     (if (< (abs (- 1 m)) (* 4 double-float-epsilon))
  86.     ;; A&S 16.6.3
  87.     (/ (lisp:cosh u))
  88.     (multiple-value-bind (v mu root-mu1)
  89.         (ascending-transform u m)
  90.       ;; A&S 16.14.4
  91.       (let* ((new-dn (elliptic-dn-ascending v mu)))
  92.         (* (/ (- 1 root-mu1) mu)
  93.            (/ (+ root-mu1 (* new-dn new-dn))
  94.           new-dn))))))
  95.  
  96.   ;; Don't use the descending version because it requires cn, dn, and
  97.   ;; sn.
  98.   (defun elliptic-cn-ascending (u m)
  99.     (if (< (abs (- 1 m)) (* 4 double-float-epsilon))
  100.     ;; A&S 16.6.2
  101.     (/ (lisp:cosh u))
  102.     (multiple-value-bind (v mu root-mu1)
  103.         (ascending-transform u m)
  104.       ;; A&S 16.14.3
  105.       (let* ((new-dn (elliptic-dn-ascending v mu)))
  106.         (* (/ (+ 1 root-mu1) mu)
  107.            (/ (- (* new-dn new-dn) root-mu1)
  108.           new-dn))))))
  109.  
  110.   ;; We don't use the ascending transform here because it requires
  111.   ;; evaluating sn, cn, and dn.  The ascending transform only needs
  112.   ;; sn.
  113.   (defun elliptic-sn-descending (u m)
  114.     ;; A&S 16.12.2
  115.     (if (< (abs m) double-float-epsilon)
  116.     (lisp:sin u)
  117.     (multiple-value-bind (v mu root-mu)
  118.         (descending-transform u m)
  119.       (let* ((new-sn (elliptic-sn-descending v mu)))
  120.         (/ (* (1+ root-mu) new-sn)
  121.            (1+ (* root-mu new-sn new-sn)))))))
  122.   #+nil
  123.   (defun elliptic-sn-ascending (u m)
  124.     (if (< (abs (- 1 m)) (* 4 double-float-epsilon))
  125.     ;; A&S 16.6.1
  126.     (tanh u)
  127.     (multiple-value-bind (v mu root-mu1)
  128.         (ascending-transform u m)
  129.       ;; A&S 16.14.2
  130.       (let* ((new-cn (elliptic-cn-ascending v mu))
  131.          (new-dn (elliptic-dn-ascending v mu))
  132.          (new-sn (elliptic-sn-ascending v mu)))
  133.         (/ (* (+ 1 root-mu1) new-sn new-cn)
  134.            new-dn)))))
  135. )
  136.  
  137. (defun sn (u m)
  138.   (if (and (realp u) (realp m))
  139.       (realpart (elliptic-sn-descending u m))
  140.       (elliptic-sn-descending u m)))
  141.  
  142. (defun cn (u m)
  143.   (if (and (realp u) (realp m))
  144.       (realpart (elliptic-cn-ascending u m))
  145.       (elliptic-cn-ascending u m)))
  146.  
  147. (defun dn (u m)
  148.   (if (and (realp u) (realp m))
  149.       (realpart (elliptic-dn-ascending u m))
  150.       (elliptic-dn-ascending u m)))
  151.  
  152.  
  153. ;;
  154. ;; How this works, I think.
  155. ;;
  156. ;; $jacobi_sn is the user visible function JACOBI_SN.  We put
  157. ;; properties on this symbol so maxima can figure out what to do with
  158. ;; it.
  159.  
  160. ;; Tell maxima how to simplify the functions $jacobi_sn, etc.  This
  161. ;; borrows heavily from trigi.lisp.
  162. (defprop %jacobi_sn simp-%jacobi_sn operators)
  163. (defprop %jacobi_cn simp-%jacobi_cn operators)
  164. (defprop %jacobi_dn simp-%jacobi_dn operators)
  165. (defprop %inverse_jacobi_sn simp-%inverse_jacobi_sn operators)
  166. (defprop %inverse_jacobi_cn simp-%inverse_jacobi_cn operators)
  167. (defprop %inverse_jacobi_dn simp-%inverse_jacobi_dn operators)
  168.  
  169. ;; Tell maxima what the derivatives are.
  170. ;;
  171. ;; Lawden says the derivative wrt to k but that's not what we want.
  172. ;;
  173. ;; Here's the derivation we used, based on how Lawden get's his results.
  174. ;;
  175. ;; Let
  176. ;;
  177. ;; diff(sn(u,m),m) = s
  178. ;; diff(cn(u,m),m) = p
  179. ;; diff(dn(u,m),m) = q
  180. ;;
  181. ;; From the derivatives of sn, cn, dn wrt to u, we have
  182. ;;
  183. ;; diff(sn(u,m),u) = cn(u)*dn(u)
  184. ;; diff(cn(u,m),u) = -cn(u)*dn(u)
  185. ;; diff(dn(u,m),u) = -m*sn(u)*cn(u)
  186. ;;
  187.  
  188. ;; Differentiate these wrt to m:
  189. ;;
  190. ;; diff(s,u) = p*dn + cn*q
  191. ;; diff(p,u) = -p*dn - q*dn
  192. ;; diff(q,u) = -sn*cn - m*s*cn - m*sn*q
  193. ;;
  194. ;; Also recall that
  195. ;;
  196. ;; sn(u)^2 + cn(u)^2 = 1
  197. ;; dn(u)^2 + m*sn(u)^2 = 1
  198. ;;
  199. ;; Differentiate these wrt to m:
  200. ;;
  201. ;; sn*s + cn*p = 0
  202. ;; 2*dn*q + sn^2 + 2*m*sn*s = 0
  203. ;;
  204. ;; Thus,
  205. ;;
  206. ;; p = -s*sn/cn
  207. ;; q = -m*s*sn/dn - sn^2/dn/2
  208. ;;
  209. ;; So
  210. ;; diff(s,u) = -s*sn*dn/cn - m*s*sn*cn/dn - sn^2*cn/dn/2
  211. ;;
  212. ;; or
  213. ;;
  214. ;; diff(s,u) + s*(sn*dn/cn + m*sn*cn/dn) = -1/2*sn^2*cn/dn
  215. ;;
  216. ;; diff(s,u) + s*sn/cn/dn*(dn^2 + m*cn^2) = -1/2*sn^2*cn/dn
  217. ;;
  218. ;; Multiply through by the integrating factor 1/cn/dn:
  219. ;;
  220. ;; diff(s/cn/dn, u) = -1/2*sn^2/dn^2 = -1/2*sd^2.
  221. ;;
  222. ;; Interate this to get
  223. ;;
  224. ;; s/cn/dn = C + -1/2*int sd^2
  225. ;;
  226. ;; It can be shown that C is zero.
  227. ;;
  228. ;; We know that (by differentiating this expression)
  229. ;;
  230. ;; int dn^2 = (1-m)*u+m*sn*cd + m*(1-m)*int sd^2
  231. ;;
  232. ;; or
  233. ;;
  234. ;; int sd^2 = 1/m/(1-m)*int dn^2 - u/m - sn*cd/(1-m)
  235. ;;
  236. ;; Thus, we get
  237. ;;
  238. ;; s/cn/dn = u/(2*m) + sn*cd/(2*(1-m)) - 1/2/m/(1-m)*int dn^2
  239. ;;
  240. ;; or
  241. ;;
  242. ;; s = 1/(2*m)*u*cn*dn + 1/(2*(1-m))*sn*cn^2 - 1/2/(m*(1-m))*cn*dn*E(u)
  243. ;;
  244. ;; where E(u) = int dn^2 = elliptic_e(am(u)) = elliptic_e(asin(sn(u)))
  245. ;;
  246. ;; This is our desired result:
  247. ;;
  248. ;; s = 1/(2*m)*cn*dn*[u - elliptic_e(asin(sn(u)),m)/(1-m)] + sn*cn^2/2/(1-m)
  249. ;;
  250. ;;
  251. ;; Since diff(cn(u,m),m) = p = -s*sn/cn, we have
  252. ;;
  253. ;; p = -1/(2*m)*sn*dn[u - elliptic_e(asin(sn(u)),m)/(1-m)] - sn^2*cn/2/(1-m)
  254. ;;
  255. ;; diff(dn(u,m),m) = q = -m*s*sn/dn - sn^2/dn/2
  256. ;;
  257. ;; q = -1/2*sn*cn*[u-elliptic_e(asin(sn),m)/(1-m)] - m*sn^2*cn^2/dn/2/(1-m)
  258. ;;
  259. ;;      - sn^2/dn/2
  260. ;;
  261. ;;   = -1/2*sn*cn*[u-elliptic_e(asin(sn),m)/(1-m)] + dn*sn^2/2/(m-1)
  262. ;;
  263. (defprop %jacobi_sn
  264.     ((u m)
  265.      ((mtimes) ((%jacobi_cn) u m) ((%jacobi_dn) u m))
  266.      ((mplus simp)
  267.       ((mtimes simp) ((rat simp) 1 2)
  268.        ((mexpt simp) ((mplus simp) 1 ((mtimes simp) -1 m)) -1)
  269.        ((mexpt simp) ((%jacobi_cn simp) u m) 2) ((%jacobi_sn simp) u m))
  270.       ((mtimes simp) ((rat simp) 1 2) ((mexpt simp) m -1)
  271.        ((%jacobi_cn simp) u m) ((%jacobi_dn simp) u m)
  272.        ((mplus simp) u
  273.     ((mtimes simp) -1 ((mexpt simp) ((mplus simp) 1 ((mtimes simp) -1 m)) -1)
  274.      (($elliptic_e simp) ((%asin simp) ((%jacobi_sn simp) u m)) m))))))
  275.   grad)
  276.  
  277. (defprop %jacobi_cn
  278.     ((u m)
  279.      ((mtimes simp) -1 ((%jacobi_sn simp) u m) ((%jacobi_dn simp) u m))
  280.      ((mplus simp)
  281.       ((mtimes simp) ((rat simp) -1 2)
  282.        ((mexpt simp) ((mplus simp) 1 ((mtimes simp) -1 m)) -1)
  283.        ((%jacobi_cn simp) u m) ((mexpt simp) ((%jacobi_sn simp) u m) 2))
  284.       ((mtimes simp) ((rat simp) -1 2) ((mexpt simp) m -1)
  285.        ((%jacobi_dn simp) u m) ((%jacobi_sn simp) u m)
  286.        ((mplus simp) u
  287.     ((mtimes simp) -1 ((mexpt simp) ((mplus simp) 1 ((mtimes simp) -1 m)) -1)
  288.      (($elliptic_e simp) ((%asin simp) ((%jacobi_sn simp) u m)) m))))))
  289.   grad)
  290.  
  291. (defprop %jacobi_dn
  292.     ((u m)
  293.      ((mtimes) -1 m ((%jacobi_sn) u m) ((%jacobi_cn) u m))
  294.      ((mplus simp)
  295.       ((mtimes simp) ((rat simp) -1 2)
  296.        ((mexpt simp) ((mplus simp) 1 ((mtimes simp) -1 m)) -1)
  297.        ((%jacobi_dn simp) u m) ((mexpt simp) ((%jacobi_sn simp) u m) 2))
  298.       ((mtimes simp) ((rat simp) -1 2) ((%jacobi_cn simp) u m)
  299.        ((%jacobi_sn simp) u m)
  300.        ((mplus simp) u
  301.     ((mtimes simp) -1
  302.      ((mexpt simp) ((mplus simp) 1 ((mtimes simp) -1 m)) -1)
  303.      (($elliptic_e simp) ((%asin simp) ((%jacobi_sn simp) u m)) m))))))
  304.   grad)
  305.  
  306. ;; The inverse elliptic functions.
  307. ;;
  308. ;; F(phi|m) = asn(sin(phi),m)
  309. ;; 
  310. ;; so asn(u,m) = F(asin(u)|m)
  311. (defprop %inverse_jacobi_sn
  312.     ((x m)
  313.      ;; 1/sqrt(1-x^2)/sqrt(1-m*x^2)
  314.      ((mtimes simp)
  315.       ((mexpt simp) ((mplus simp) 1 ((mtimes simp) -1 ((mexpt simp) x 2)))
  316.        ((rat simp) -1 2))
  317.       ((mexpt simp) ((mplus simp) 1 ((mtimes simp) -1 m ((mexpt simp) x 2)))
  318.        ((rat simp) -1 2)))
  319.      ;; diff(F(asin(u)|m),m)
  320.      ((mtimes simp) ((mexpt simp) ((mplus simp) 1 ((mtimes simp) -1 m)) -1)
  321.       ((mplus simp)
  322.        ((mtimes simp) -1 x
  323.     ((mexpt simp) ((mplus simp) 1 ((mtimes simp) -1 ((mexpt simp) x 2)))
  324.      ((rat simp) 1 2))
  325.     ((mexpt simp) ((mplus simp) 1 ((mtimes simp) -1 m ((mexpt simp) x 2)))
  326.      ((rat simp) -1 2)))
  327.        ((mtimes simp) ((mexpt simp) m -1)
  328.     ((mplus simp) ((%elliptic_e simp) ((%asin simp) x) m)
  329.      ((mtimes simp) -1 ((mplus simp) 1 ((mtimes simp) -1 m))
  330.       ((%elliptic_f simp) ((%asin simp) x) m)))))))
  331.   grad)
  332.  
  333. ;; Let u = inverse_jacobi_cn(x).  Then jacobi_cn(u) = x or
  334. ;; sqrt(1-jacobi_sn(u)^2) = x.  Or
  335. ;;
  336. ;; jacobi_sn(u) = sqrt(1-x^2)
  337. ;;
  338. ;; So u = inverse_jacobi_sn(sqrt(1-x^2),m) = inverse_jacob_cn(x,m)
  339. ;;
  340. (defprop %inverse_jacobi_cn
  341.     ((x m)
  342.      ;; -1/sqrt(1-x^2)/sqrt(1-m*x^2)
  343.      ((mtimes simp) -1
  344.       ((mexpt simp) ((mplus simp) 1 ((mtimes simp) -1 ((mexpt simp) x 2)))
  345.        ((rat simp) -1 2))
  346.       ((mexpt simp) ((mplus simp) 1 ((mtimes simp) -1 m ((mexpt simp) x 2)))
  347.        ((rat simp) -1 2)))
  348.      ((mtimes simp) ((mexpt simp) ((mplus simp) 1 ((mtimes simp) -1 m)) -1)
  349.       ((mplus simp)
  350.        ((mtimes simp) -1
  351.     ((mexpt simp)
  352.      ((mplus simp) 1
  353.       ((mtimes simp) -1 m ((mplus simp) 1 ((mtimes simp) -1 ((mexpt simp) x 2)))))
  354.      ((rat simp) -1 2))
  355.     ((mexpt simp) ((mplus simp) 1 ((mtimes simp) -1 ((mexpt simp) x 2))) ((rat simp) 1 2))
  356.     ((mabs simp) x))
  357.        ((mtimes simp) ((mexpt simp) m -1)
  358.     ((mplus simp)
  359.      ((%elliptic_e simp)
  360.       ((%asin simp)
  361.        ((mexpt simp) ((mplus simp) 1 ((mtimes simp) -1 ((mexpt simp) x 2))) ((rat simp) 1 2)))
  362.       m)
  363.      ((mtimes simp) -1 ((mplus simp) 1 ((mtimes simp) -1 m))
  364.       ((%elliptic_f simp)
  365.        ((%asin simp)
  366.         ((mexpt simp) ((mplus simp) 1 ((mtimes simp) -1 ((mexpt simp) x 2))) ((rat simp) 1 2)))
  367.        m)))))))
  368.   grad)
  369.  
  370. ;; Let u = inverse_jacobi_dn(x).  Then
  371. ;;
  372. ;; jacobi_dn(u) = x or
  373. ;;
  374. ;; x^2 = jacobi_dn(u)^2 = 1 - m*jacobi_sn(u)^2
  375. ;;
  376. ;; so jacobi_sn(u) = sqrt(1-x^2)/sqrt(m)
  377. ;;
  378. ;; or u = inverse_jacobi_sn(sqrt(1-x^2)/sqrt(m))
  379. (defprop %inverse_jacobi_dn
  380.     ((x m)
  381.      ;; -1/sqrt(1-x^2)/sqrt(x^2+m-1)
  382.      ((mtimes simp)
  383.       ((mexpt simp) ((mplus simp) 1 ((mtimes simp) -1 ((mexpt simp) x 2)))
  384.        ((rat simp) -1 2))
  385.       ((mexpt simp) ((mplus simp) -1 m ((mexpt simp) x 2)) ((rat simp) -1 2)))
  386.      ((mplus simp)
  387.       ((mtimes simp) ((rat simp) -1 2) ((mexpt simp) m ((rat simp) -3 2))
  388.        ((mexpt simp)
  389.     ((mplus simp) 1
  390.      ((mtimes simp) -1 ((mexpt simp) m -1)
  391.       ((mplus simp) 1 ((mtimes simp) -1 ((mexpt simp) x 2)))))
  392.     ((rat simp) -1 2))
  393.        ((mexpt simp) ((mplus simp) 1 ((mtimes simp) -1 ((mexpt simp) x 2)))
  394.     ((rat simp) 1 2))
  395.        ((mexpt simp) ((mabs simp) x) -1))
  396.       ((mtimes simp) ((mexpt simp) ((mplus simp) 1 ((mtimes simp) -1 m)) -1)
  397.        ((mplus simp)
  398.     ((mtimes simp) -1 ((mexpt simp) m ((rat simp) -1 2))
  399.      ((mexpt simp)
  400.       ((mplus simp) 1
  401.        ((mtimes simp) -1 ((mexpt simp) m -1)
  402.         ((mplus simp) 1 ((mtimes simp) -1 ((mexpt simp) x 2)))))
  403.       ((rat simp) 1 2))
  404.      ((mexpt simp) ((mplus simp) 1 ((mtimes simp) -1 ((mexpt simp) x 2)))
  405.       ((rat simp) 1 2))
  406.      ((mexpt simp) ((mabs simp) x) -1))
  407.     ((mtimes simp) ((mexpt simp) m -1)
  408.      ((mplus simp)
  409.       ((%elliptic_e simp)
  410.        ((%asin simp)
  411.         ((mtimes simp) ((mexpt simp) m ((rat simp) -1 2))
  412.          ((mexpt simp) ((mplus simp) 1
  413.                 ((mtimes simp) -1 ((mexpt simp) x 2)))
  414.           ((rat simp) 1 2))))
  415.        m)
  416.       ((mtimes simp) -1 ((mplus simp) 1 ((mtimes simp) -1 m))
  417.        ((%elliptic_f simp)
  418.         ((%asin simp)
  419.          ((mtimes simp) ((mexpt simp) m ((rat simp) -1 2))
  420.           ((mexpt simp) ((mplus simp) 1
  421.                  ((mtimes simp) -1 ((mexpt simp) x 2)))
  422.            ((rat simp) 1 2))))
  423.         m))))))))
  424.   grad)
  425.  
  426. ;; Define the actual functions for the user
  427. (defmfun $jacobi_sn (u m)
  428.   (simplify (list '(%jacobi_sn) (resimplify u) (resimplify m))))
  429. (defmfun $jacobi_cn (u m)
  430.   (simplify (list '(%jacobi_cn) (resimplify u) (resimplify m))))
  431. (defmfun $jacobi_dn (u m)
  432.   (simplify (list '(%jacobi_dn) (resimplify u) (resimplify m))))
  433.  
  434. (defmfun $inverse_jacobi_sn (u m)
  435.   (simplify (list '(%inverse_jacobi_sn) (resimplify u) (resimplify m))))
  436.  
  437. (defmfun $inverse_jacobi_cn (u m)
  438.   (simplify (list '(%inverse_jacobi_cn) (resimplify u) (resimplify m))))
  439.  
  440. (defmfun $inverse_jacobi_dn (u m)
  441.   (simplify (list '(%inverse_jacobi_dn) (resimplify u) (resimplify m))))
  442.  
  443. ;; A complex number looks like
  444. ;;
  445. ;; ((MPLUS SIMP) 0.70710678118654757 ((MTIMES SIMP) 0.70710678118654757 $%I))
  446. ;;
  447. ;; or
  448. ;;
  449. ;; ((MPLUS SIMP) 1 $%I))
  450. ;;
  451. (defun complex-number-p (u)
  452.   ;; Return non-NIL if U is a complex number (or number)
  453.   (or (numberp u)
  454.       (and (consp u)
  455.        (numberp (second u))
  456.        (or (and (consp (third u))
  457.             (numberp (second (third u)))
  458.             (eq (third (third u)) '$%i))
  459.            (and (eq (third u) '$%i))))))
  460.  
  461. (defun complexify (x)
  462.   ;; Convert a Lisp number to a maxima number
  463.   (if (realp x)
  464.       x
  465.       (if (zerop (realpart x))
  466.       `((mtimes) ,(imagpart x) $%i)
  467.       `((mplus simp) ,(realpart x)
  468.         ((mtimes simp) ,(imagpart x) $%i)))))
  469.  
  470.  
  471. (defun kc-arg (exp m)
  472.   ;; Replace elliptic_kc(m) in the expression with sym.  Check to see
  473.   ;; if the resulting expression is linear in sym and the constant
  474.   ;; term is zero.  If so, return the coefficient of sym, i.e, the
  475.   ;; coefficient of elliptic_kc(m).
  476.   (let* ((sym (gensym))
  477.      (arg (maxima-substitute sym `((%elliptic_kc) ,m) exp)))
  478.     (if (and (not (equalp arg exp))
  479.          (linearp arg sym)
  480.          (zerop1 (coefficient arg sym 0)))
  481.     (coefficient arg sym 1)
  482.     nil)))
  483.  
  484. (defun kc-arg2 (exp m)
  485.   ;; Replace elliptic_kc(m) in the expression with sym.  Check to see
  486.   ;; if the resulting expression is linear in sym and the constant
  487.   ;; term is zero.  If so, return the coefficient of sym, i.e, the
  488.   ;; coefficient of elliptic_kc(m), and the constant term.  Otherwise,
  489.   ;; return NIL.
  490.   (let* ((sym (gensym))
  491.      (arg (maxima-substitute sym `((%elliptic_kc) ,m) exp)))
  492.     (if (and (not (equalp arg exp))
  493.          (linearp arg sym))
  494.     (list (coefficient arg sym 1)
  495.           (coefficient arg sym 0))
  496.     nil)))
  497.  
  498. ;; Tell maxima how to simplify the functions
  499. ;;
  500. ;; FORM is list containing the actual expression.  I don't really know
  501. ;; what Y and Z contain.  Most of this modeled after SIMP-%SIN.
  502. (defmfun simp-%jacobi_sn (form y z)
  503.   (declare (ignore y))
  504.   (twoargcheck form)
  505.   (let ((u (simpcheck (cadr form) z))
  506.     (m (simpcheck (caddr form) z))
  507.     coef)
  508.     (cond ((or (and (floatp u) (floatp m))
  509.            (and $numer (numberp u) (numberp m)))
  510.        ;; Numerically evaluate sn
  511.        (sn (float u 1d0) (float m 1d0)))
  512.       ((and $numer (complex-number-p u)
  513.         (complex-number-p m))
  514.        ;; For complex values.  Should we really do this?
  515.        (let ((result (sn (complex ($realpart u) ($imagpart u))
  516.                  (complex ($realpart m) ($imagpart m)))))
  517.          (complexify result)))
  518.       ((zerop1 u)
  519.        ;; A&S 16.5.1
  520.        0)
  521.       ((zerop1 m)
  522.        ;; A&S 16.6.1
  523.        `((%sin) ,u))
  524.       ((onep1 m)
  525.        ;; A&S 16.6.1
  526.        `((%tanh) ,u))
  527.       ((and $trigsign (mminusp* u))
  528.        (neg (cons-exp '%jacobi_sn (neg u) m)))
  529.       ;; A&S 16.20.1 (Jacobi's Imaginary transformation)
  530.       ((and $%iargs (multiplep u '$%i))
  531.        (mul '$%i
  532.         (cons-exp '%jacobi_sc (coeff u '$%i 1) (add 1 (neg m)))))
  533.       ((setq coef (kc-arg2 u m))
  534.        ;; sn(m*K+u)
  535.        ;;
  536.        ;; A&S 16.8.1
  537.        (destructuring-bind (lin const)
  538.            coef
  539.          (cond ((integerp lin)
  540.             (ecase (mod lin 4)
  541.               (0
  542.                ;; sn(4*m*K + u) = sn(u), sn(0) = 0
  543.                (if (zerop1 const)
  544.                0
  545.                `((%jacobi_sn simp) ,const ,m)))
  546.               (1
  547.                ;; sn(4*m*K + K + u) = sn(K+u) = cd(u)
  548.                ;; sn(K) = 1
  549.                (if (zerop1 const)
  550.                1
  551.                `((%jacobi_cd simp) ,const ,m)))
  552.               (2
  553.                ;; sn(4*m*K+2*K + u) = sn(2*K+u) = -sn(u)
  554.                ;; sn(2*K) = 0
  555.                (if (zerop1 const)
  556.                0
  557.                (neg `((%jacobi_sn simp) ,const ,m))))
  558.               (3
  559.                ;; sn(4*m*K+3*K+u) = sn(2*K + K + u) = -sn(K+u) = -cd(u)
  560.                ;; sn(3*K) = -1
  561.                (if (zerop1 const)
  562.                -1
  563.                (neg `((%jacobi_cd simp) ,const ,m))))))
  564.            ((and (alike1 lin 1//2)
  565.              (zerop1 const))
  566.             ;; A&S 16.5.2
  567.             ;;
  568.             ;; sn(1/2*K) = 1/sqrt(1+sqrt(1-m))
  569.             `((mexpt simp)
  570.               ((mplus simp) 1
  571.                ((mexpt simp)
  572.             ((mplus simp) 1 ((mtimes simp) -1 ,m))
  573.             ((rat simp) 1 2)))
  574.               ((rat) -1 2)))
  575.            ((and (alike1 lin 3//2)
  576.              (zerop1 const))
  577.             ;; A&S 16.5.2
  578.             ;;
  579.             ;; sn(1/2*K + K) = cd(1/2*K,m)
  580.             (simplifya
  581.              `((%jacobi_cd) ((mtimes) ((rat) 1 2) ((%elliptic_kc) ,m))
  582.                ,m)
  583.              nil))
  584.            (t
  585.             (eqtest (list '(%jacobi_sn) u m) form)))))
  586.       (t
  587.        ;; Nothing to do
  588.        (eqtest (list '(%jacobi_sn) u m) form)))))
  589.  
  590. (defmfun simp-%jacobi_cn (form y z)
  591.   (declare (ignore y))
  592.   (twoargcheck form)
  593.   (let ((u (simpcheck (cadr form) z))
  594.     (m (simpcheck (caddr form) z))
  595.     coef)
  596.     (cond ((or (and (floatp u) (floatp m))
  597.            (and $numer (numberp u) (numberp m)))
  598.        (cn (float u 1d0) (float m 1d0)))
  599.       ((and $numer (complex-number-p u)
  600.         (complex-number-p m))
  601.        (let ((result (cn (complex ($realpart u) ($imagpart u))
  602.                  (complex ($realpart m) ($imagpart m)))))
  603.          (complexify result)))
  604.       ((zerop1 u)
  605.        ;; A&S 16.5.1
  606.        1)
  607.       ((zerop1 m)
  608.        ;; A&S 16.6.2
  609.        `((%cos) ,u))
  610.       ((onep1 m)
  611.        ;; A&S 16.6.2
  612.        `((%sech) ,u))
  613.       ((and $trigsign (mminusp* u))
  614.        (cons-exp '%jacobi_cn (neg u) m))
  615.       ;; A&S 16.20.2 (Jacobi's Imaginary transformation)
  616.       ((and $%iargs (multiplep u '$%i))
  617.        (cons-exp '%jacobi_nc (coeff u '$%i 1) (add 1 (neg m))))
  618.       ((setq coef (kc-arg2 u m))
  619.        ;; cn(m*K+u)
  620.        ;;
  621.        ;; A&S 16.8.2
  622.        (destructuring-bind (lin const)
  623.            coef
  624.          (cond ((integerp lin)
  625.             (ecase (mod lin 4)
  626.               (0
  627.                ;; cn(4*m*K + u) = cn(u),
  628.                ;; cn(0) = 1
  629.                (if (zerop1 const)
  630.                1
  631.                `((%jacobi_cn simp) ,const ,m)))
  632.               (1
  633.                ;; cn(4*m*K + K + u) = cn(K+u) = -sqrt(m1)*sd(u)
  634.                ;; cn(K) = 0
  635.                (if (zerop1 const)
  636.                0
  637.                (neg `((mtimes simp)
  638.                   ((mexpt simp)
  639.                    ((mplus simp) 1 ((mtimes simp) -1 ,m))
  640.                    ((rat simp) 1 2))
  641.                   ((%jacobi_sd simp) ,const ,m)))))
  642.               (2
  643.                ;; cn(4*m*K + 2*K + u) = cn(2*K+u) = -cn(u)
  644.                ;; cn(2*K) = -1
  645.                (if (zerop1 const)
  646.                -1
  647.                (neg `((%jacobi_cn) ,const ,m))))
  648.               (3
  649.                ;; cn(4*m*K + 3*K + u) = cn(2*K + K + u) =
  650.                ;; -cn(K+u) = sqrt(m1)*sd(u)
  651.                ;;
  652.                ;; cn(3*K) = 0
  653.                (if (zerop1 const)
  654.                0
  655.                `((mtimes simp)
  656.                  ((mexpt simp)
  657.                   ((mplus simp) 1 ((mtimes simp) -1 ,m))
  658.                   ((rat simp) 1 2))
  659.                  ((%jacobi_sd simp) ,const ,m))))))
  660.            ((and (alike1 lin 1//2)
  661.              (zerop1 const))
  662.             ;; A&S 16.5.2
  663.             ;; cn(1/2*K) = (1-m)^(1/4)/sqrt(1+sqrt(1-m))
  664.             `((mtimes simp)
  665.               ((mexpt simp) ((mplus simp) 1
  666.                      ((mtimes simp) -1 ,m))
  667.                ((rat simp) 1 4))
  668.               ((mexpt simp)
  669.                ((mplus simp) 1
  670.             ((mexpt simp)
  671.              ((mplus simp) 1
  672.               ((mtimes simp) -1 ,m))
  673.              ((rat simp) 1 2)))
  674.                ((rat simp) -1 2))))
  675.            (t
  676.             (eqtest (list '(%jacobi_cn) u m) form)))))
  677.       (t
  678.        (eqtest (list '(%jacobi_cn) u m) form)))))
  679.  
  680. (defmfun simp-%jacobi_dn (form y z)
  681.   (declare (ignore y))
  682.   (twoargcheck form)
  683.   (let ((u (simpcheck (cadr form) z))
  684.     (m (simpcheck (caddr form) z))
  685.     coef)
  686.     (cond ((or (and (floatp u) (floatp m))
  687.            (and $numer (numberp u) (numberp m)))
  688.        (dn (float u 1d0) (float m 1d0)))
  689.       ((and $numer (complex-number-p u)
  690.         (complex-number-p m))
  691.        (let ((result (dn (complex ($realpart u) ($imagpart u))
  692.                  (complex ($realpart m) ($imagpart m)))))
  693.          (complexify result)))
  694.       ((zerop1 u)
  695.        ;; A&S 16.5.1
  696.        1)
  697.       ((zerop1 m)
  698.        ;; A&S 16.6.3
  699.        1)
  700.       ((onep1 m)
  701.        ;; A&S 16.6.3
  702.        `(($sech) ,u))
  703.       ((and $trigsign (mminusp* u))
  704.        (cons-exp '%jacobi_dn (neg u) m))
  705.       ;; A&S 16.20.2 (Jacobi's Imaginary transformation)
  706.       ((and $%iargs (multiplep u '$%i))
  707.        (cons-exp '%jacobi_dc (coeff u '$%i 1)
  708.              (add 1 (neg m))))
  709.       ((setq coef (kc-arg2 u m))
  710.        ;; A&S 16.8.3
  711.        ;;
  712.        ;; dn(m*K+u) has period 2K
  713.        ;;
  714.        (destructuring-bind (lin const)
  715.            coef
  716.          (cond ((integerp lin)
  717.             (ecase (mod lin 2)
  718.               (0
  719.                ;; dn(2*m*K + u) = dn(u)
  720.                ;; dn(0) = 1
  721.                (if (zerop1 const)
  722.                1
  723.                ;; dn(4*m*K+2*K + u) = dn(2*K+u) = dn(u)
  724.                `((%jacobi_dn) ,const ,m)))
  725.               (1
  726.                ;; dn(2*m*K + K + u) = dn(K + u) = sqrt(1-m)*nd(u)
  727.                ;; dn(K) = sqrt(1-m)
  728.                (if (zerop1 const)
  729.                `((mexpt simp)
  730.                  ((mplus simp) 1 ((mtimes simp) -1 ,m))
  731.                  ((rat simp) 1 2))
  732.                `((mtimes simp)
  733.                  ((mexpt simp)
  734.                   ((mplus simp) 1 ((mtimes simp) -1 ,m))
  735.                   ((rat simp) 1 2))
  736.                  ((%jacobi_nd simp) ,const ,m))))))
  737.            ((and (alike1 lin 1//2)
  738.              (zerop1 const))
  739.             ;; A&S 16.5.2
  740.             ;; dn(1/2*K) = (1-m)^(1/4)
  741.             `((mexpt simp)
  742.               ((mplus simp) 1 ((mtimes simp) -1 ,m))
  743.               ((rat simp) 1 4)))
  744.            (t
  745.             (eqtest (list '(%jacobi_cn) u m) form)))))
  746.       (t (eqtest (list '(%jacobi_dn) u m) form)))))
  747.  
  748. ;; Should we simplify the inverse elliptic functions into the
  749. ;; appropriate incomplete elliptic integral?  I think we should leave
  750. ;; it, but perhaps allow some way to do that transformation if
  751. ;; desired.
  752.  
  753. (defmfun simp-%inverse_jacobi_sn (form y z)
  754.   (declare (ignore y))
  755.   (twoargcheck form)
  756.   (let ((u (simpcheck (cadr form) z))
  757.     (m (simpcheck (caddr form) z)))
  758.     (cond ((or (and (floatp u) (floatp m))
  759.            (and $numer (numberp u) (numberp m)))
  760.        ;; Numerically evaluate asn
  761.        ;;
  762.        ;; asn(x,m) = F(asin(x),m)
  763.        (elliptic-f (lisp:asin u) m))
  764.       ((zerop1 u)
  765.        ;; asn(0,m) = 0
  766.        0)
  767.       ((onep1 u)
  768.        ;; asn(1,m) = elliptic_kc(m)
  769.        `(($elliptic_kc) ,m))
  770.       ((and (numberp u) (onep1 (- u)))
  771.        ;; asn(-1,m) = -elliptic_kc(m)
  772.        `((mtimes) -1 ((%elliptic_kc) ,m)))
  773.       ((zerop1 m)
  774.        ;; asn(x,0) = F(asin(x),0) = asin(x)
  775.        `((%elliptic_f) ((%asin) ,u) 0))
  776.       ((onep1 m)
  777.        ;; asn(x,1) = F(asin(x),1) = log(tan(pi/2+asin(x)/2))
  778.        `((%elliptic_f) ((%asin) ,u) 1))
  779.       (t
  780.        ;; Nothing to do
  781.        (eqtest (list '(%inverse_jacobi_sn) u m) form)))))
  782.  
  783. (defmfun simp-%inverse_jacobi_cn (form y z)
  784.   (declare (ignore y))
  785.   (twoargcheck form)
  786.   (let ((u (simpcheck (cadr form) z))
  787.     (m (simpcheck (caddr form) z)))
  788.     (cond ((or (and (floatp u) (floatp m))
  789.            (and $numer (numberp u) (numberp m)))
  790.        ;; Numerically evaluate acn
  791.        ;;
  792.        ;; acn(x,m) = F(acos(x),m)
  793.        (elliptic-f (acos u) m))
  794.       ((zerop1 m)
  795.        ;; asn(x,0) = F(acos(x),0) = acos(x)
  796.        `((%elliptic_f) ((%acos) ,u) 0))
  797.       ((onep1 m)
  798.        ;; asn(x,1) = F(asin(x),1) = log(tan(pi/2+asin(x)/2))
  799.        `((%elliptic_f) ((%acos) ,u) 1))
  800.       ((zerop1 u)
  801.        `((%elliptic_kc) ,m))
  802.       ((onep1 u)
  803.        0)
  804.       (t
  805.        ;; Nothing to do
  806.        (eqtest (list '(%inverse_jacobi_cn) u m) form)))))
  807.  
  808. ;;
  809. ;; adn(u) = x.
  810. ;; u = dn(x) = sqrt(1-m*sn(x)^2)
  811. ;; sn(x)^2 = (1-u^2)/m
  812. ;; sn(x) = sqrt((1-u^2)/m)
  813. ;; x = asn(sqrt((1-u^2)/m))
  814. ;; x = adn(u) = asn(sqrt((1-u^2)/m))
  815.  
  816. (defmfun simp-%inverse_jacobi_dn (form y z)
  817.   (declare (ignore y))
  818.   (twoargcheck form)
  819.   (let ((u (simpcheck (cadr form) z))
  820.     (m (simpcheck (caddr form) z)))
  821.     (cond ((or (and (floatp u) (floatp m))
  822.            (and $numer (numberp u) (numberp m)))
  823.        ;; Numerically evaluate adn
  824.        (let ((phi (/ (* (sqrt (- 1 u)) (sqrt (+ 1 u)))
  825.              (sqrt m))))
  826.          (elliptic-f (asin phi) m)))
  827.       ((onep1 m)
  828.        ;; x = dn(u,1) = sech(u).  so u = asech(x)
  829.        `((%asech) ,u))
  830.       (t
  831.        ;; Nothing to do
  832.        (eqtest (list '(%inverse_jacobi_dn) u m) form)))))
  833.  
  834. ;;;; Elliptic integrals
  835.  
  836. ;; Carlson's elliptic integral of the first kind.
  837. ;;
  838. ;; ***PURPOSE  Compute the incomplete or complete elliptic integral of the
  839. ;;             1st kind.  For X, Y, and Z non-negative and at most one of
  840. ;;             them zero, RF(X,Y,Z) = Integral from zero to infinity of
  841. ;;                                 -1/2     -1/2     -1/2
  842. ;;                       (1/2)(t+X)    (t+Y)    (t+Z)    dt.
  843. ;;             If X, Y or Z is zero, the integral is complete.
  844. ;;
  845. ;;          Value of IER assigned by the DRF routine
  846. ;;
  847. ;;                   Value assigned         Error Message Printed
  848. ;;                   IER = 1                MIN(X,Y,Z) .LT. 0.0D0
  849. ;;                       = 2                MIN(X+Y,X+Z,Y+Z) .LT. LOLIM
  850. ;;                       = 3                MAX(X,Y,Z) .GT. UPLIM
  851. ;;
  852. ;;    Special double precision functions via DRF
  853. ;;
  854. ;;
  855. ;;
  856. ;;
  857. ;;                   Legendre form of ELLIPTIC INTEGRAL of 1st kind
  858. ;;
  859. ;;                   -----------------------------------------
  860. ;;
  861. ;;
  862. ;;
  863. ;;                                              2         2   2
  864. ;;                   F(PHI,K) = SIN(PHI) DRF(COS (PHI),1-K SIN (PHI),1)
  865. ;;
  866. ;;
  867. ;;                                   2
  868. ;;                   K(K) = DRF(0,1-K ,1)
  869. ;;
  870. ;;
  871. ;;                          PI/2     2   2      -1/2
  872. ;;                        = INT  (1-K SIN (PHI) )   D PHI
  873. ;;                           0
  874. ;;
  875. ;;
  876. ;;
  877. ;;                   Bulirsch form of ELLIPTIC INTEGRAL of 1st kind
  878. ;;
  879. ;;                   -----------------------------------------
  880. ;;
  881. ;;
  882. ;;                                           2 2    2
  883. ;;                   EL1(X,KC) = X DRF(1,1+KC X ,1+X )
  884. ;;
  885. ;;
  886. ;;                   Lemniscate constant A
  887. ;;
  888. ;;                   -----------------------------------------
  889. ;;
  890. ;;
  891. ;;                        1      4 -1/2
  892. ;;                   A = INT (1-S )    DS = DRF(0,1,2) = DRF(0,2,1)
  893. ;;                        0
  894. ;;
  895. ;;
  896. ;;
  897. ;;     -------------------------------------------------------------------
  898. ;;
  899. ;; ***REFERENCES  B. C. Carlson and E. M. Notis, Algorithms for incomplete
  900. ;;                  elliptic integrals, ACM Transactions on Mathematical
  901. ;;                  Software 7, 3 (September 1981), pp. 398-403.
  902. ;;                B. C. Carlson, Computing elliptic integrals by
  903. ;;                  duplication, Numerische Mathematik 33, (1979),
  904. ;;                  pp. 1-16.
  905. ;;                B. C. Carlson, Elliptic integrals of the first kind,
  906. ;;                  SIAM Journal of Mathematical Analysis 8, (1977),
  907. ;;                 pp. 231-242.
  908.  
  909. (let ((errtol (expt (* 4 double-float-epsilon) 1/6))
  910.       (uplim (/ most-positive-double-float 5))
  911.       (lolim (* #-gcl least-positive-normalized-double-float
  912.         #+gcl least-positive-double-float
  913.         5))
  914.       (c1 (float 1/24 1d0))
  915.       (c2 (float 3/44 1d0))
  916.       (c3 (float 1/14 1d0)))
  917.   (declare (double-float errtol c1 c2 c3))
  918.   (defun drf (x y z)
  919.     "Compute Carlson's incomplete or complete elliptic integral of the
  920. first kind:
  921.  
  922.                    INF
  923.                   /
  924.                   [                     1
  925.   RF(x, y, z) =   I    ----------------------------------- dt
  926.                   ]    SQRT(x + t) SQRT(y + t) SQRT(z + t)
  927.                   /
  928.                    0
  929.  
  930.  
  931. where x >= 0, y >= 0, z >=0, and at most one of x, y, z is zero.
  932. "
  933.     (declare (double-float x y z))
  934.     ;; Check validity of input
  935.     (assert (and (>= x 0) (>= y 0) (>= z 0)
  936.          (plusp (+ x y)) (plusp (+ x z)) (plusp (+ y z))))
  937.     (assert (<= (max x y z) uplim))
  938.     (assert (>= (min (+ x y) (+ x z) (+ y z)) lolim))
  939.     (locally 
  940.     (declare (type (double-float 0d0) x y)
  941.          (type (double-float (0d0)) z)
  942.          (optimize (speed 3)))
  943.       (loop
  944.       (let* ((mu (/ (+ x y z) 3))
  945.          (x-dev (- 2 (/ (+ mu x) mu)))
  946.          (y-dev (- 2 (/ (+ mu y) mu)))
  947.          (z-dev (- 2 (/ (+ mu z) mu))))
  948.         (when (< (max (abs x-dev) (abs y-dev) (abs z-dev)) errtol)
  949.           (let ((e2 (- (* x-dev y-dev) (* z-dev z-dev)))
  950.             (e3 (* x-dev y-dev z-dev)))
  951.         (return (/ (+ 1
  952.                   (* e2 (- (* c1 e2)
  953.                        1/10
  954.                        (* c2 e3)))
  955.                   (* c3 e3))
  956.                (sqrt mu)))))
  957.         (let* ((x-root (sqrt x))
  958.            (y-root (sqrt y))
  959.            (z-root (sqrt z))
  960.            (lam (+ (* x-root (+ y-root z-root)) (* y-root z-root))))
  961.           (setf x (* (+ x lam) 1/4))
  962.           (setf y (* (+ y lam) 1/4))
  963.           (setf z (* (+ z lam) 1/4))))))))
  964.  
  965. ;; Elliptic integral of the first kind (Legendre's form):
  966. ;;
  967. ;;
  968. ;;      phi
  969. ;;     /
  970. ;;     [             1
  971. ;;     I    ------------------- ds
  972. ;;     ]                  2
  973. ;;     /    SQRT(1 - m SIN (s))
  974. ;;     0
  975.  
  976. (defun elliptic-f (phi-arg m-arg)
  977.   (let ((phi (float phi-arg 1d0))
  978.     (m (float m-arg 1d0)))
  979.     (cond ((> m 1)
  980.        ;; A&S 17.4.15
  981.        (/ (elliptic-f (asin (* (sqrt m) (sin phi))) (/ m))))
  982.       ((< m 0)
  983.        ;; A&S 17.4.17
  984.        (let* ((m (- m))
  985.           (m+1 (+ 1 m))
  986.           (root (sqrt m+1))
  987.           (m/m+1 (/ m m+1)))
  988.          (- (/ (elliptic-f (float (/ pi 2) 1d0) m/m+1)
  989.            root)
  990.         (/ (elliptic-f (- (float (/ pi 2) 1d0) phi) m/m+1)
  991.            root))))
  992.       ((= m 0)
  993.        ;; A&S 17.4.19
  994.        phi)
  995.       ((= m 1)
  996.        ;; A&S 17.4.21
  997.        (log (lisp:tan (+ (/ phi 2) (float (/ pi 2) 1d0)))))
  998.       ((minusp phi)
  999.        (- (elliptic-f (- phi) m)))
  1000.       ((> phi pi)
  1001.        ;; A&S 17.4.3
  1002.        (multiple-value-bind (s phi-rem)
  1003.            (truncate phi (float pi 1d0))
  1004.          (+ (* 2 s (elliptic-k m))
  1005.         (elliptic-f phi-rem m))))
  1006.       ((<= phi (/ pi 2))
  1007.        (let ((sin-phi (sin phi))
  1008.          (cos-phi (cos phi))
  1009.          (k (sqrt m)))
  1010.          (* sin-phi
  1011.         (drf (* cos-phi cos-phi)
  1012.              (* (- 1 (* k sin-phi))
  1013.             (+ 1 (* k sin-phi)))
  1014.              1d0))))
  1015.       ((< phi pi)
  1016.        (+ (* 2 (elliptic-k m))
  1017.           (elliptic-f (- phi (float pi 1d0)) m))))))
  1018.  
  1019. ;; Complete elliptic integral of the first kind
  1020. (defun elliptic-k (m)
  1021.   (declare (double-float m))
  1022.   (cond ((> m 1)
  1023.      ;; A&S 17.4.15
  1024.      (/ (elliptic-f (asin (sqrt m)) (/ m))))
  1025.     ((< m 0)
  1026.      ;; A&S 17.4.17
  1027.      (let* ((m (- m))
  1028.         (m+1 (+ 1 m))
  1029.         (root (sqrt m+1))
  1030.         (m/m+1 (/ m m+1)))
  1031.        (- (/ (elliptic-k m/m+1)
  1032.          root)
  1033.           (/ (elliptic-f 0d0 m/m+1)
  1034.          root))))
  1035.     ((= m 0)
  1036.      ;; A&S 17.4.19
  1037.      (float (/ pi 2) 1d0))
  1038.     (t
  1039.      (let ((k (sqrt m)))
  1040.        (drf 0d0 (* (- 1 k)
  1041.                (+ 1 k))
  1042.         1d0)))))
  1043. ;;
  1044. ;; Carlsons' elliptic integral of the second kind.
  1045. ;;
  1046. ;;   1.     DRD
  1047. ;;          Evaluate an INCOMPLETE (or COMPLETE) ELLIPTIC INTEGRAL
  1048. ;;          of the second kind
  1049. ;;          Standard FORTRAN function routine
  1050. ;;          Double precision version
  1051. ;;          The routine calculates an approximation result to
  1052. ;;          DRD(X,Y,Z) = Integral from zero to infinity of
  1053. ;;                              -1/2     -1/2     -3/2
  1054. ;;                    (3/2)(t+X)    (t+Y)    (t+Z)    dt,
  1055. ;;          where X and Y are nonnegative, X + Y is positive, and Z is
  1056. ;;          positive.  If X or Y is zero, the integral is COMPLETE.
  1057.  
  1058. ;;    -------------------------------------------------------------------
  1059. ;;
  1060. ;;
  1061. ;;   Special double precision functions via DRD and DRF
  1062. ;;
  1063. ;;
  1064. ;;                  Legendre form of ELLIPTIC INTEGRAL of 2nd kind
  1065. ;;
  1066. ;;                  -----------------------------------------
  1067. ;;
  1068. ;;
  1069. ;;                                             2         2   2
  1070. ;;                  E(PHI,K) = SIN(PHI) DRF(COS (PHI),1-K SIN (PHI),1) -
  1071. ;;
  1072. ;;                     2      3             2         2   2
  1073. ;;                  -(K/3) SIN (PHI) DRD(COS (PHI),1-K SIN (PHI),1)
  1074. ;;
  1075. ;;
  1076. ;;                                  2        2            2
  1077. ;;                  E(K) = DRF(0,1-K ,1) - (K/3) DRD(0,1-K ,1)
  1078. ;;
  1079. ;;                         PI/2     2   2      1/2
  1080. ;;                       = INT  (1-K SIN (PHI) )  D PHI
  1081. ;;                          0
  1082. ;;
  1083. ;;                  Bulirsch form of ELLIPTIC INTEGRAL of 2nd kind
  1084. ;;
  1085. ;;                  -----------------------------------------
  1086. ;;
  1087. ;;                                               2 2    2
  1088. ;;                  EL2(X,KC,A,B) = AX DRF(1,1+KC X ,1+X ) +
  1089. ;;
  1090. ;;                                              3          2 2    2
  1091. ;;                                 +(1/3)(B-A) X DRD(1,1+KC X ,1+X )
  1092. ;;
  1093. ;;
  1094. ;;
  1095. ;;
  1096. ;;                  Legendre form of alternative ELLIPTIC INTEGRAL
  1097. ;;                  of 2nd kind
  1098. ;;
  1099. ;;                  -----------------------------------------
  1100. ;;
  1101. ;;
  1102. ;;
  1103. ;;                            Q     2       2   2  -1/2
  1104. ;;                  D(Q,K) = INT SIN P  (1-K SIN P)     DP
  1105. ;;                            0
  1106. ;;
  1107. ;;
  1108. ;;
  1109. ;;                                     3          2     2   2
  1110. ;;                  D(Q,K) = (1/3) (SIN Q) DRD(COS Q,1-K SIN Q,1)
  1111. ;;
  1112. ;;
  1113. ;;
  1114. ;;
  1115. ;;                  Lemniscate constant  B
  1116. ;;
  1117. ;;                  -----------------------------------------
  1118. ;;
  1119. ;;
  1120. ;;
  1121. ;;
  1122. ;;                       1    2    4 -1/2
  1123. ;;                  B = INT  S (1-S )    DS
  1124. ;;                       0
  1125. ;;
  1126. ;;
  1127. ;;                  B = (1/3) DRD (0,2,1)
  1128. ;;
  1129. ;;
  1130. ;;                  Heuman's LAMBDA function
  1131. ;;
  1132. ;;                  -----------------------------------------
  1133. ;;
  1134. ;;
  1135. ;;
  1136. ;;                  (PI/2) LAMBDA0(A,B) =
  1137. ;;
  1138. ;;                                    2                2
  1139. ;;                 = SIN(B) (DRF(0,COS (A),1)-(1/3) SIN (A) *
  1140. ;;
  1141. ;;                            2               2         2       2
  1142. ;;                  *DRD(0,COS (A),1)) DRF(COS (B),1-COS (A) SIN (B),1)
  1143. ;;
  1144. ;;                            2       3             2
  1145. ;;                  -(1/3) COS (A) SIN (B) DRF(0,COS (A),1) *
  1146. ;;
  1147. ;;                           2         2       2
  1148. ;;                   *DRD(COS (B),1-COS (A) SIN (B),1)
  1149. ;;
  1150. ;;
  1151. ;;
  1152. ;;                  Jacobi ZETA function
  1153. ;;
  1154. ;;                  -----------------------------------------
  1155. ;;
  1156. ;;                             2                 2       2   2
  1157. ;;                  Z(B,K) = (K/3) SIN(B) DRF(COS (B),1-K SIN (B),1)
  1158. ;;
  1159. ;;
  1160. ;;                                       2             2
  1161. ;;                             *DRD(0,1-K ,1)/DRF(0,1-K ,1)
  1162. ;;
  1163. ;;                               2       3           2       2   2
  1164. ;;                            -(K /3) SIN (B) DRD(COS (B),1-K SIN (B),1)
  1165. ;;
  1166. ;;
  1167.  
  1168. (let ((errtol (expt (/ double-float-epsilon 3) 1/6))
  1169.       (c1 (float 3/14 1d0))
  1170.       (c2 (float 1/6 1d0))
  1171.       (c3 (float 9/22 1d0))
  1172.       (c4 (float 3/26 1d0)))
  1173.   (declare (double-float errtol c1 c2 c3 c4))
  1174.   (defun drd (x y z)
  1175.     (declare (real x y z))
  1176.     ;; Check validity of input
  1177.  
  1178.     (assert (and (>= x 0) (>= y 0) (>= z 0) (plusp (+ x y)))
  1179.         (x y z))
  1180.         
  1181.     (let ((x (float x 1d0))
  1182.       (y (float y 1d0))
  1183.       (z (float z 1d0))
  1184.       (sigma 0d0)
  1185.       (power4 1d0))
  1186.       (declare (type (double-float 0d0) x y power4 sigma)
  1187.            (type (double-float (0d0)) z)
  1188.            (optimize (speed 3)))
  1189.       (loop
  1190.       (let* ((mu (* 1/5 (+ x y (* 3 z))))
  1191.          (x-dev (/ (- mu x) mu))
  1192.          (y-dev (/ (- mu y) mu))
  1193.          (z-dev (/ (- mu z) mu)))
  1194.         (when (< (max (abs x-dev) (abs y-dev) (abs z-dev)) errtol)
  1195.           (let* ((ea (* x-dev y-dev))
  1196.              (eb (* z-dev z-dev))
  1197.              (ec (- ea eb))
  1198.              (ed (- ea (* 6 eb)))
  1199.              (ef (+ ed ec ec))
  1200.              (s1 (* ed (+ (- c1)
  1201.                   (* 1/4 c3 ed)
  1202.                   (* -3/2 c4 z-dev ef))))
  1203.              (s2 (* z-dev (+ (* c2 ef)
  1204.                      (* z-dev (+ (* (- c3) ec)
  1205.                          (* z-dev c4 ea)))))))
  1206.         (return (+ (* 3 sigma)
  1207.                (/ (* power4 (+ 1 s1 s2))
  1208.                   (* mu (sqrt mu)))))))
  1209.         (let* ((x-root (sqrt x))
  1210.            (y-root (sqrt y))
  1211.            (z-root (sqrt z))
  1212.            (lam (+ (* x-root (+ y-root z-root)) (* y-root z-root))))
  1213.           (incf sigma (/ power4 z-root (+ z lam)))
  1214.           (setf power4 (/ power4 4))
  1215.           (setf x (/ (+ x lam) 4))
  1216.           (setf y (/ (+ y lam) 4))
  1217.           (setf z (/ (+ z lam) 4))))))))
  1218.  
  1219. ;; Elliptic integral of the second kind (Legendre's form):
  1220. ;;
  1221. ;;
  1222. ;;      phi
  1223. ;;     /
  1224. ;;     [                  2
  1225. ;;     I    SQRT(1 - m SIN (s)) ds
  1226. ;;     ]
  1227. ;;     /
  1228. ;;      0
  1229.  
  1230. (defun elliptic-e (phi m)
  1231.   (declare (double-float phi m))
  1232.   (cond ((= m 0)
  1233.      ;; A&S 17.4.23
  1234.      phi)
  1235.     ((= m 1)
  1236.      ;; A&S 17.4.25
  1237.      (sin phi))
  1238.     (t
  1239.      (let* ((sin-phi (sin phi))
  1240.         (cos-phi (cos phi))
  1241.         (k (sqrt m))
  1242.         (y (* (- 1 (* k sin-phi))
  1243.               (+ 1 (* k sin-phi)))))
  1244.        (- (* sin-phi
  1245.          (drf (* cos-phi cos-phi) y 1d0))
  1246.           (* (/ m 3)
  1247.          (expt sin-phi 3)
  1248.          (drd (* cos-phi cos-phi) y 1d0)))))))
  1249.  
  1250. ;; Complete version
  1251. (defun elliptic-ec (m)
  1252.   (declare (double-float m))
  1253.   (cond ((= m 0)
  1254.      ;; A&S 17.4.23
  1255.      (float (/ pi 2) 1d0))
  1256.     ((= m 1)
  1257.      ;; A&S 17.4.25
  1258.      1d0)
  1259.     (t
  1260.      (let* ((k (sqrt m))
  1261.         (y (* (- 1 k)
  1262.               (+ 1 k))))
  1263.        (- (drf 0d0 y 1d0)
  1264.           (* (/ m 3)
  1265.          (drd 0d0 y 1d0)))))))
  1266.  
  1267.  
  1268. ;; Define the elliptic integrals for maxima
  1269. ;;
  1270. ;; We use the definitions given in A&S 17.2.6 and 17.2.8.  In particular:
  1271. ;;
  1272. ;;                 phi
  1273. ;;                /
  1274. ;;                [             1
  1275. ;; F(phi|m)  =    I    ------------------- ds
  1276. ;;                ]                  2
  1277. ;;                /    SQRT(1 - m SIN (s))
  1278. ;;                 0
  1279. ;;
  1280. ;; and
  1281. ;;
  1282. ;;              phi
  1283. ;;             /
  1284. ;;             [                  2
  1285. ;; E(phi|m) =  I    SQRT(1 - m SIN (s)) ds
  1286. ;;             ]
  1287. ;;             /
  1288. ;;              0
  1289. ;;
  1290. ;; That is, we do not use the modular angle, alpha, as the second arg;
  1291. ;; the parameter m = sin(alpha)^2 is used.
  1292. ;;
  1293.  
  1294.  
  1295. (defprop $elliptic_f simp-$elliptic_f operators)
  1296. (defprop $elliptic_e simp-$elliptic_e operators)
  1297.  
  1298. ;; The derivative of F(phi|m) wrt to phi is easy.  The derivative wrt
  1299. ;; to m is harder.  Here is a derivation.  Hope I got it right.
  1300. ;;
  1301. ;; diff(integrate(1/sqrt(1-m*sin(x)^2),x,0,phi), m);
  1302. ;;
  1303. ;;                PHI
  1304. ;;               /           2
  1305. ;;               [        SIN (x)
  1306. ;;               I    ------------------ dx
  1307. ;;               ]         2    3/2
  1308. ;;               /    (1 - m SIN (x))
  1309. ;;                0
  1310. ;;                --------------------------
  1311. ;;                       2
  1312. ;;
  1313. ;; 
  1314. ;; Now use the following relationship that is easily verified:
  1315. ;;
  1316. ;;               2                 2
  1317. ;;    (1 - m) SIN (x)           COS (x)                 COS(x) SIN(x)
  1318. ;;  ------------------- = ------------------- - DIFF(-------------------, x)
  1319. ;;                 2                     2                          2
  1320. ;;  SQRT(1 - m SIN (x))   SQRT(1 - m SIN (x))         SQRT(1 - m SIN (x))
  1321. ;;
  1322. ;;
  1323. ;; Now integrate this to get:
  1324. ;;
  1325. ;; 
  1326. ;;            PHI
  1327. ;;           /            2
  1328. ;;           [         SIN (x)
  1329. ;;    (1 - m) I       ------------------- dx =
  1330. ;;           ]             2
  1331. ;;           /       SQRT(1 - m SIN (x))
  1332. ;;            0
  1333.  
  1334. ;;
  1335. ;;                PHI
  1336. ;;               /            2
  1337. ;;               [         COS (x)
  1338. ;;              + I    ------------------- dx
  1339. ;;               ]             2
  1340. ;;               /    SQRT(1 - m SIN (x))
  1341. ;;                0
  1342. ;;                    COS(PHI) SIN(PHI)
  1343. ;;                -  ---------------------
  1344. ;;                        2
  1345. ;;                  SQRT(1 - m SIN (PHI))
  1346. ;;
  1347. ;; Use the fact that cos(x)^2 = 1 - sin(x)^2 to show that this
  1348. ;; integral on the RHS is:
  1349. ;;
  1350. ;;
  1351. ;;          (1 - m) elliptic_F(PHI, m) + elliptic_E(PHI, m)
  1352. ;;           -------------------------------------------
  1353. ;;                       m
  1354. ;; So, finally, we have
  1355. ;;
  1356. ;;
  1357. ;; 
  1358. ;;   d                
  1359. ;; 2 -- (elliptic_F(PHI, m)) = 
  1360. ;;   dm                
  1361. ;;
  1362. ;;  elliptic_E(PHI, m) - (1 - m) elliptic_F(PHI, m)     COS(PHI) SIN(PHI)
  1363. ;;  ---------------------------------------------- - ---------------------
  1364. ;;                m                      2
  1365. ;;                              SQRT(1 - m SIN (PHI))
  1366. ;;   ----------------------------------------------------------------------
  1367. ;;                      1 - m
  1368.  
  1369. (defprop $elliptic_f
  1370.     ((phi m)
  1371.      ;; diff wrt phi
  1372.      ;; 1/sqrt(1-m*sin(phi)^2)
  1373.      ((mexpt simp)
  1374.       ((mplus simp) 1 ((mtimes simp) -1 m ((mexpt simp) ((%sin simp) phi) 2)))
  1375.       ((rat simp) -1 2))
  1376.      ;; diff wrt m
  1377.      ((mtimes simp) ((rat simp) 1 2)
  1378.       ((mexpt simp) ((mplus simp) 1 ((mtimes simp) -1 m)) -1)
  1379.       ((mplus simp)
  1380.        ((mtimes simp) ((mexpt simp) m -1)
  1381.     ((mplus simp) (($elliptic_e simp) phi m)
  1382.      ((mtimes simp) -1 ((mplus simp) 1 ((mtimes simp) -1 m))
  1383.       (($elliptic_f simp) phi m))))
  1384.        ((mtimes simp) -1 ((%cos simp) phi) ((%sin simp) phi)
  1385.     ((mexpt simp)
  1386.      ((mplus simp) 1
  1387.       ((mtimes simp) -1 m ((mexpt simp) ((%sin simp) phi) 2)))
  1388.      ((rat simp) -1 2))))))
  1389.   grad)
  1390.  
  1391. ;;
  1392. ;; The derivative of E(phi|m) wrt to m is much simpler to derive than F(phi|m).
  1393. ;;
  1394. ;; Take the derivative of the definition to get
  1395. ;;
  1396. ;;         PHI
  1397. ;;        /         2
  1398. ;;        [          SIN (x)
  1399. ;;        I    ------------------- dx
  1400. ;;        ]              2
  1401. ;;        /    SQRT(1 - m SIN (x))
  1402. ;;         0
  1403. ;;      - ---------------------------
  1404. ;;             2
  1405. ;;
  1406. ;; It is easy to see that
  1407. ;;
  1408. ;;                 PHI
  1409. ;;                /         2
  1410. ;;                [          SIN (x)
  1411. ;;  elliptic_F(PHI, m) - m I    ------------------- dx = elliptic_E(PHI, m)
  1412. ;;                ]              2
  1413. ;;                /    SQRT(1 - m SIN (x))
  1414. ;;                 0
  1415. ;;
  1416. ;; So we finally have
  1417. ;;
  1418. ;;   d                   elliptic_E(PHI, m) - elliptic_F(PHI, m)
  1419. ;;   -- (elliptic_E(PHI, m)) = ---------------------------------------
  1420. ;;   dm                            2 m
  1421.  
  1422. (defprop $elliptic_e
  1423.     ((phi m)
  1424.      ;; (1-m*sin(phi)^2)
  1425.      ((mplus simp) 1 ((mtimes simp) -1 m ((mexpt simp) ((%sin simp) phi) 2)))
  1426.      ;; diff wrt m
  1427.      ((mtimes simp) ((rat simp) 1 2) ((mexpt simp) m -1)
  1428.       ((mplus simp) (($elliptic_e simp) phi m)
  1429.        ((mtimes simp) -1 (($elliptic_f simp) phi m)))))
  1430.   grad)
  1431.             
  1432. (defmfun $elliptic_f (phi m)
  1433.   (simplify (list '($elliptic_f) (resimplify phi) (resimplify m))))
  1434. (defmfun $elliptic_e (phi m)
  1435.   (simplify (list '($elliptic_e) (resimplify phi) (resimplify m))))
  1436.  
  1437. (defmfun simp-$elliptic_f (form y z)
  1438.   (declare (ignore y))
  1439.   (twoargcheck form)
  1440.   (let ((phi (simpcheck (cadr form) z))
  1441.     (m (simpcheck (caddr form) z)))
  1442.     (cond ((or (and (floatp phi) (floatp m))
  1443.            (and $numer (numberp phi) (numberp m)))
  1444.        ;; Numerically evaluate it
  1445.        (elliptic-f (float phi 1d0) (float m 1d0)))
  1446.       ((zerop1 phi)
  1447.        0)
  1448.       ((zerop1 m)
  1449.        ;; A&S 17.4.19
  1450.        phi)
  1451.       ((onep1 m)
  1452.        ;; A&S 17.4.21.  Let's pick the log tan form.
  1453.        `((%log) ((%tan)
  1454.              ((mplus) ((mtimes) $%pi ((rat) 1 4))
  1455.               ((mtimes) ((rat) 1 2) ,phi)))))
  1456.       ((alike1 phi '((mtimes) ((rat) 1 2) $%pi))
  1457.        ;; Complete elliptic integral
  1458.        `((%elliptic_kc) ,m))
  1459.       (t
  1460.        ;; Nothing to do
  1461.        (eqtest (list '($elliptic_f) phi m) form)))))
  1462.  
  1463. (defmfun simp-$elliptic_e (form y z)
  1464.   (declare (ignore y))
  1465.   (twoargcheck form)
  1466.   (let ((phi (simpcheck (cadr form) z))
  1467.     (m (simpcheck (caddr form) z)))
  1468.     (cond ((or (and (floatp phi) (floatp m))
  1469.            (and $numer (numberp phi) (numberp m)))
  1470.        ;; Numerically evaluate it
  1471.        (elliptic-e (float phi 1d0) (float m 1d0)))
  1472.       ((zerop1 phi)
  1473.        0)
  1474.       ((zerop1 m)
  1475.        ;; A&S 17.4.23
  1476.        phi)
  1477.       ((onep1 m)
  1478.        ;; A&S 17.4.25
  1479.        `((%sin) ,phi))
  1480.       ((alike1 phi '((mtimes) ((rat) 1 2) $%pi))
  1481.        ;; Complete elliptic integral
  1482.        `((%elliptic_ec) ,m))
  1483.       (t
  1484.        ;; Nothing to do
  1485.        (eqtest (list '($elliptic_e) phi m) form)))))
  1486.     
  1487.  
  1488. ;; Complete elliptic integrals
  1489. ;;
  1490. ;; elliptic_kc(m) = elliptic_f(%pi/2, m)
  1491. ;;
  1492. ;; elliptic_ec(m) = elliptic_e(%pi/2, m)
  1493. ;;
  1494. (defmfun $elliptic_kc (m)
  1495.   (simplify (list '(%elliptic_kc) (resimplify m))))
  1496. (defmfun $elliptic_ec (m)
  1497.   (simplify (list '(%elliptic_ec) (resimplify m))))
  1498.  
  1499.  
  1500. (defprop %elliptic_kc simp-%elliptic_kc operators)
  1501. (defprop %elliptic_ec simp-%elliptic_ec operators)
  1502.  
  1503. (defmfun simp-%elliptic_kc (form y z)
  1504.   (declare (ignore y))
  1505.   (oneargcheck form)
  1506.   (let ((m (simpcheck (cadr form) z)))
  1507.     (cond ((or (and (floatp m))
  1508.            (and $numer (numberp m)))
  1509.        ;; Numerically evaluate it
  1510.        (elliptic-k (float m 1d0)))
  1511.       ((zerop1 m)
  1512.        '((mtimes) ((rat) 1 2) $%pi))
  1513.       (t
  1514.        ;; Nothing to do
  1515.        (eqtest (list '(%elliptic_kc) m) form)))))
  1516.  
  1517. (defprop %elliptic_kc
  1518.     ((m)
  1519.      ;; diff wrt m
  1520.      ((mtimes)
  1521.       ((mplus) ((%elliptic_ec) m)
  1522.        ((mtimes) -1
  1523.     ((%elliptic_kc) m)
  1524.     ((mplus) 1 ((mtimes) -1 m))))
  1525.       ((mexpt) ((mplus) 1 ((mtimes) -1 m)) -1)
  1526.       ((mexpt) m -1)))
  1527.   grad)
  1528.  
  1529. (defmfun simp-%elliptic_ec (form y z)
  1530.   (declare (ignore y))
  1531.   (oneargcheck form)
  1532.   (let ((m (simpcheck (cadr form) z)))
  1533.     (cond ((or (and (floatp m))
  1534.            (and $numer (numberp m)))
  1535.        ;; Numerically evaluate it
  1536.        (elliptic-ec (float m 1d0)))
  1537.       ((zerop1 m)
  1538.        '((mtimes) ((rat) 1 2) $%pi))
  1539.       ;; Some special cases we know about.
  1540.       (t
  1541.        ;; Nothing to do
  1542.        (eqtest (list '(%elliptic_ec) m) form)))))
  1543.  
  1544. (defprop %elliptic_ec
  1545.     ((m)
  1546.      ((mtimes) ((rat) 1 2)
  1547.       ((mplus) ((%elliptic_ec) m)
  1548.        ((mtimes) -1 ((%elliptic_kc)
  1549.                m)))
  1550.       ((mexpt) m -1)))
  1551.   grad)
  1552.  
  1553. ;;
  1554. ;; Elliptic integral of the third kind:
  1555. ;;
  1556. ;; (A&S 17.2.14)
  1557. ;;
  1558. ;;                 phi
  1559. ;;                /
  1560. ;;                [                     1
  1561. ;; PI(n;phi|m) =  I    ----------------------------------- ds
  1562. ;;                ]                  2               2
  1563. ;;                /    SQRT(1 - m SIN (s)) (1 - n SIN (s))
  1564. ;;                 0
  1565. ;;
  1566. ;; As with E and F, we do not use the modular angle alpha but the
  1567. ;; parameter m = sin(alpha)^2.
  1568. ;;
  1569. (defprop $elliptic_pi simp-$elliptic_pi operators)
  1570.  
  1571. (defmfun $elliptic_pi (n phi m)
  1572.   (simplify (list '($elliptic_pi)
  1573.           (resimplify n) (resimplify phi) (resimplify m))))
  1574.  
  1575. (defmfun simp-$elliptic_pi (form y z)
  1576.   (declare (ignore y))
  1577.   ;;(threeargcheck form)
  1578.   (let ((n (simpcheck (cadr form) z))
  1579.     (phi (simpcheck (caddr form) z))
  1580.     (m (simpcheck (cadddr form) z)))
  1581.     (cond ((or (and (floatp n) (floatp phi) (floatp m))
  1582.            (and $numer (numberp n) (numberp phi) (numberp m)))
  1583.        ;; Numerically evaluate it
  1584.        (elliptic-pi (float n 1d0) (float phi 1d0) (float m 1d0)))
  1585.       ((zerop1 n)
  1586.        `(($elliptic_f) ,phi ,m))
  1587.       ((zerop1 m)
  1588.        ;; 3 cases depending on n < 1, n > 1, or n = 1.
  1589.        (let ((s (asksign `((mplus) -1 ,n))))
  1590.          (case s
  1591.            ($positive
  1592.         `((mtimes)
  1593.           ((mexpt) ((mplus) -1 ,n) ((rat) -1 2))
  1594.           ((%atanh)
  1595.            ((mtimes) ((mexpt) ((mplus) -1 ,n) ((rat) 1 2))
  1596.             ((%tan) ,phi)))))
  1597.            ($negative
  1598.         `((mtimes)
  1599.           ((mexpt) ((mplus) 1 ((mtimes) -1 ,n)) ((rat) -1 2))
  1600.           ((%atan)
  1601.            ((mtimes) ((mexpt)
  1602.                   ((mplus) 1 ((mtimes) -1 ,n))
  1603.                   ((rat) 1 2))
  1604.             ((%tan) ,phi)))))
  1605.            ($zero
  1606.         `((%tan) ,phi)))))
  1607.       (t
  1608.        ;; Nothing to do
  1609.        (eqtest (list '($elliptic_pi) n phi m) form)))))
  1610.  
  1611. (defun elliptic-pi (n phi m)
  1612.   ;; Note: Carlson's DRJ has n defined as the negative of the n given
  1613.   ;; in A&S.
  1614.   (let* ((nn (- n))
  1615.     (sin-phi (sin phi))
  1616.     (cos-phi (cos phi))
  1617.     (k (sqrt m))
  1618.     (k2sin (* (- 1 (* k sin-phi))
  1619.           (+ 1 (* k sin-phi)))))
  1620.     (- (* sin-phi (drf (expt cos-phi 2) k2sin 1d0))
  1621.        (* (/ nn 3) (expt sin-phi 3)
  1622.       (drj (expt cos-phi 2) k2sin 1d0
  1623.            (- 1 (* n (expt sin-phi 2))))))))
  1624.     
  1625. ;;***PURPOSE  Calculate a double precision approximation to
  1626. ;;             DRC(X,Y) = Integral from zero to infinity of
  1627. ;;                              -1/2     -1
  1628. ;;                    (1/2)(t+X)    (t+Y)  dt,
  1629. ;;            where X is nonnegative and Y is positive.
  1630. ;;   1.     DRC
  1631. ;;          Standard FORTRAN function routine
  1632. ;;          Double precision version
  1633. ;;          The routine calculates an approximation result to
  1634. ;;          DRC(X,Y) = integral from zero to infinity of
  1635. ;;
  1636. ;;                              -1/2     -1
  1637. ;;                    (1/2)(t+X)    (t+Y)  dt,
  1638. ;;
  1639. ;;          where X is nonnegative and Y is positive.  The duplication
  1640. ;;          theorem is iterated until the variables are nearly equal,
  1641. ;;          and the function is then expanded in Taylor series to fifth
  1642. ;;          order.  Logarithmic, inverse circular, and inverse hyper-
  1643. ;;          bolic functions can be expressed in terms of DRC.
  1644. ;;
  1645. ;;   --------------------------------------------------------------------
  1646. ;;
  1647. ;;   Special functions via DRC
  1648. ;;
  1649. ;;
  1650. ;;
  1651. ;;                  LN X                X .GT. 0
  1652. ;;
  1653. ;;                                             2
  1654. ;;                  LN(X) = (X-1) DRC(((1+X)/2)  , X )
  1655. ;;
  1656. ;;
  1657. ;;   --------------------------------------------------------------------
  1658. ;;
  1659. ;;                  ARCSIN X            -1 .LE. X .LE. 1
  1660. ;;
  1661. ;;                                       2
  1662. ;;                  ARCSIN X = X DRC (1-X  ,1 )
  1663. ;;
  1664. ;;   --------------------------------------------------------------------
  1665. ;;
  1666. ;;                  ARCCOS X            0 .LE. X .LE. 1
  1667. ;;
  1668. ;;
  1669. ;;                                     2       2
  1670. ;;                  ARCCOS X = SQRT(1-X ) DRC(X  ,1 )
  1671. ;;
  1672. ;;   --------------------------------------------------------------------
  1673. ;;
  1674. ;;                  ARCTAN X            -INF .LT. X .LT. +INF
  1675. ;;
  1676. ;;                                        2
  1677. ;;                  ARCTAN X = X DRC(1,1+X  )
  1678. ;;
  1679. ;;   --------------------------------------------------------------------
  1680. ;;
  1681. ;;                  ARCCOT X            0 .LE. X .LT. INF
  1682. ;;
  1683. ;;                                  2   2
  1684. ;;                  ARCCOT X = DRC(X  ,X +1 )
  1685. ;;
  1686. ;;   --------------------------------------------------------------------
  1687. ;;
  1688. ;;                  ARCSINH X           -INF .LT. X .LT. +INF
  1689. ;;
  1690. ;;                                       2
  1691. ;;                  ARCSINH X = X DRC(1+X  ,1 )
  1692. ;;
  1693. ;;   --------------------------------------------------------------------
  1694. ;;
  1695. ;;                  ARCCOSH X           X .GE. 1
  1696. ;;
  1697. ;;                                    2         2
  1698. ;;                  ARCCOSH X = SQRT(X -1) DRC(X  ,1 )
  1699. ;;
  1700. ;;   --------------------------------------------------------------------
  1701. ;;
  1702. ;;                  ARCTANH X           -1 .LT. X .LT. 1
  1703. ;;
  1704. ;;                                         2
  1705. ;;                  ARCTANH X = X DRC(1,1-X  )
  1706. ;;
  1707. ;;   --------------------------------------------------------------------
  1708. ;;
  1709. ;;                  ARCCOTH X           X .GT. 1
  1710. ;;
  1711. ;;                                   2   2
  1712. ;;                  ARCCOTH X = DRC(X  ,X -1 )
  1713. ;;
  1714. ;;   --------------------------------------------------------------------
  1715. ;;
  1716. ;;***REFERENCES  B. C. Carlson and E. M. Notis, Algorithms for incomplete
  1717. ;;                 elliptic integrals, ACM Transactions on Mathematical
  1718. ;;                 Software 7, 3 (September 1981), pp. 398-403.
  1719. ;;               B. C. Carlson, Computing elliptic integrals by
  1720. ;;                 duplication, Numerische Mathematik 33, (1979),
  1721. ;;                 pp. 1-16.
  1722. ;;               B. C. Carlson, Elliptic integrals of the first kind,
  1723. ;;                 SIAM Journal of Mathematical Analysis 8, (1977),
  1724. ;;                 pp. 231-242.
  1725.  
  1726.  
  1727. (let ((errtol (expt (/ double-float-epsilon 3) 1/6))
  1728.       (c1 (float 1/7 1d0))
  1729.       (c2 (float 9/22 1d0)))
  1730.   (declare (double-float errtol c1 c2))
  1731.   (defun drc (x y)
  1732.     (declare (type (double-float 0d0) x)
  1733.          (type (double-float (0d0)) y)
  1734.          (optimize (speed 3)))
  1735.     (let ((xn x)
  1736.       (yn y))
  1737.       (declare (type (double-float 0d0) xn)
  1738.            (type (double-float (0d0)) yn))
  1739.       (loop
  1740.       (let* ((mu (/ (+ xn yn yn) 3))
  1741.          (sn (- (/ (+ yn mu) mu) 2)))
  1742.         (declare (type double-float sn))
  1743.         (when (< (abs sn) errtol)
  1744.           (let ((s (* sn sn
  1745.               (+ 0.3d0
  1746.                  (* sn (+ c1 (* sn (+ 0.375d0 (* sn c2)))))))))
  1747.         (return-from drc (/ (+ 1 s) (sqrt mu)))))
  1748.         (let ((lam (+ (* 2 (sqrt xn) (sqrt yn)) yn)))
  1749.           (setf xn (* (+ xn lam) 0.25d0))
  1750.           (setf yn (* (+ yn lam) 0.25d0))))))))
  1751.  
  1752. ;;   1.     DRJ
  1753. ;;          Standard FORTRAN function routine
  1754. ;;          Double precision version
  1755. ;;          The routine calculates an approximation result to
  1756. ;;          DRJ(X,Y,Z,P) = Integral from zero to infinity of
  1757. ;;
  1758. ;;                                -1/2     -1/2     -1/2     -1
  1759. ;;                      (3/2)(t+X)    (t+Y)    (t+Z)    (t+P)  dt,
  1760. ;;
  1761. ;;          where X, Y, and Z are nonnegative, at most one of them is
  1762. ;;          zero, and P is positive.  If X or Y or Z is zero, the
  1763. ;;          integral is COMPLETE.  The duplication theorem is iterated
  1764. ;;          until the variables are nearly equal, and the function is
  1765. ;;          then expanded in Taylor series to fifth order.
  1766. ;;
  1767. ;;
  1768. ;;    -------------------------------------------------------------------
  1769. ;;
  1770. ;;
  1771. ;;   Special double precision functions via DRJ and DRF
  1772. ;;
  1773. ;;
  1774. ;;                  Legendre form of ELLIPTIC INTEGRAL of 3rd kind
  1775. ;;                  -----------------------------------------
  1776. ;;
  1777. ;;
  1778. ;;                          PHI         2         -1
  1779. ;;             P(PHI,K,N) = INT (1+N SIN (THETA) )   *
  1780. ;;                           0
  1781. ;;
  1782. ;;
  1783. ;;                                  2    2         -1/2
  1784. ;;                             *(1-K  SIN (THETA) )     D THETA
  1785. ;;
  1786. ;;
  1787. ;;                                           2          2   2
  1788. ;;                        = SIN (PHI) DRF(COS (PHI), 1-K SIN (PHI),1)
  1789. ;;
  1790. ;;                                   3             2         2   2
  1791. ;;                         -(N/3) SIN (PHI) DRJ(COS (PHI),1-K SIN (PHI),
  1792. ;;
  1793. ;;                                  2
  1794. ;;                         1,1+N SIN (PHI))
  1795. ;;
  1796. ;;
  1797. ;;
  1798. ;;                  Bulirsch form of ELLIPTIC INTEGRAL of 3rd kind
  1799. ;;                  -----------------------------------------
  1800. ;;
  1801. ;;
  1802. ;;                                            2 2    2
  1803. ;;                  EL3(X,KC,P) = X DRF(1,1+KC X ,1+X ) +
  1804. ;;
  1805. ;;                                            3           2 2    2     2
  1806. ;;                               +(1/3)(1-P) X  DRJ(1,1+KC X ,1+X ,1+PX )
  1807. ;;
  1808. ;;
  1809. ;;                                           2
  1810. ;;                  CEL(KC,P,A,B) = A RF(0,KC ,1) +
  1811. ;;
  1812. ;;
  1813. ;;                                                      2
  1814. ;;                                 +(1/3)(B-PA) DRJ(0,KC ,1,P)
  1815. ;;
  1816. ;;
  1817. ;;                  Heuman's LAMBDA function
  1818. ;;                  -----------------------------------------
  1819. ;;
  1820. ;;
  1821. ;;                                2                      2      2    1/2
  1822. ;;                  L(A,B,P) =(COS (A)SIN(B)COS(B)/(1-COS (A)SIN (B))   )
  1823. ;;
  1824. ;;                                            2         2       2
  1825. ;;                            *(SIN(P) DRF(COS (P),1-SIN (A) SIN (P),1)
  1826. ;;
  1827. ;;                                 2       3            2       2
  1828. ;;                            +(SIN (A) SIN (P)/(3(1-COS (A) SIN (B))))
  1829. ;;
  1830. ;;                                    2         2       2
  1831. ;;                            *DRJ(COS (P),1-SIN (A) SIN (P),1,1-
  1832. ;;
  1833. ;;                                2       2          2       2
  1834. ;;                            -SIN (A) SIN (P)/(1-COS (A) SIN (B))))
  1835. ;;
  1836. ;;
  1837. ;;
  1838. ;;                  (PI/2) LAMBDA0(A,B) =L(A,B,PI/2) =
  1839. ;;
  1840. ;;                   2                         2       2    -1/2
  1841. ;;              = COS (A)  SIN(B) COS(B) (1-COS (A) SIN (B))
  1842. ;;
  1843. ;;                           2                  2       2
  1844. ;;                 *DRF(0,COS (A),1) + (1/3) SIN (A) COS (A)
  1845. ;;
  1846. ;;                                      2       2    -3/2
  1847. ;;                 *SIN(B) COS(B) (1-COS (A) SIN (B))
  1848. ;;
  1849. ;;                           2         2       2          2       2
  1850. ;;                 *DRJ(0,COS (A),1,COS (A) COS (B)/(1-COS (A) SIN (B)))
  1851. ;;
  1852. ;;
  1853. ;;                  Jacobi ZETA function
  1854. ;;                  -----------------------------------------
  1855. ;;
  1856. ;;                        2                     2   2    1/2
  1857. ;;             Z(B,K) = (K/3) SIN(B) COS(B) (1-K SIN (B))
  1858. ;;
  1859. ;;
  1860. ;;                                  2      2   2                 2
  1861. ;;                        *DRJ(0,1-K ,1,1-K SIN (B)) / DRF (0,1-K ,1)
  1862. ;;
  1863. ;;
  1864. ;;
  1865. ;;***REFERENCES  B. C. Carlson and E. M. Notis, Algorithms for incomplete
  1866. ;;                 elliptic integrals, ACM Transactions on Mathematical
  1867. ;;                 Software 7, 3 (September 1981), pp. 398-403.
  1868. ;;               B. C. Carlson, Computing elliptic integrals by
  1869. ;;                 duplication, Numerische Mathematik 33, (1979),
  1870. ;;                 pp. 1-16.
  1871. ;;               B. C. Carlson, Elliptic integrals of the first kind,
  1872. ;;                 SIAM Journal of Mathematical Analysis 8, (1977),
  1873. ;;                 pp. 231-242.
  1874. (let ((errtol (expt (/ double-float-epsilon 3) 1/6))
  1875.       (c1 3/14)
  1876.       (c2 1/3)
  1877.       (c3 3/22)
  1878.       (c4 3/26))
  1879.   (defun drj (x y z p)
  1880.     (declare (type (double-float 0d0) x y z)
  1881.          (type (double-float (0d0)) p))
  1882.     (let ((xn x)
  1883.       (yn y)
  1884.       (zn z)
  1885.       (pn p)
  1886.       (sigma 0d0)
  1887.       (power4 1d0))
  1888.       (declare (type (double-float 0d0) xn yn zn)
  1889.            (type (double-float (0d0)) pn)
  1890.            (type double-float sigma power4))
  1891.       (loop
  1892.       (let* ((mu (* 0.2d0 (+ xn yn zn pn pn)))
  1893.          (xndev (/ (- mu xn) mu))
  1894.          (yndev (/ (- mu yn) mu))
  1895.          (zndev (/ (- mu zn) mu))
  1896.          (pndev (/ (- mu pn) mu))
  1897.          (eps (max (abs xndev) (abs yndev) (abs zndev) (abs pndev))))
  1898.         (when (< eps errtol)
  1899.           (let* ((ea (+ (* xndev (+ yndev zndev))
  1900.                 (* yndev zndev)))
  1901.              (eb (* xndev yndev zndev))
  1902.              (ec (* pndev pndev))
  1903.              (e2 (- ea (* 3 ec)))
  1904.              (e3 (+ eb (* 2 pndev (- ea ec))))
  1905.              (s1 (+ 1 (* e2 (+ (- c1)
  1906.                        (* 3/4 c3 e2)
  1907.                        (* -3/2 c4 e3)))))
  1908.              (s2 (* eb (+ (* 1/2 c2)
  1909.                   (* pndev (+ (- c3) (- c3)
  1910.                           (* pndev c4))))))
  1911.              (s3 (- (* pndev ea (- c2 (* pndev c3)))
  1912.                 (* c2 pndev ec))))
  1913.         (return-from drj (+ (* 3 sigma)
  1914.                     (/ (* power4 (+ s1 s2 s3))
  1915.                        (* mu (sqrt mu)))))))
  1916.         (let* ((xnroot (sqrt xn))
  1917.            (ynroot (sqrt yn))
  1918.            (znroot (sqrt zn))
  1919.            (lam (+ (* xnroot (+ ynroot znroot))
  1920.                (* ynroot znroot)))
  1921.            (alfa (expt (+ (* pn (+ xnroot ynroot znroot))
  1922.                   (* xnroot ynroot znroot))
  1923.                    2))
  1924.            (beta (* pn (expt (+ pn lam) 2))))
  1925.           (incf sigma (* power4 (drc alfa beta)))
  1926.           (setf power4 (* power4 0.25d0))
  1927.           (setf xn (* 0.25d0 (+ xn lam)))
  1928.           (setf yn (* 0.25d0 (+ yn lam)))
  1929.           (setf zn (* 0.25d0 (+ zn lam)))
  1930.           (setf pn (* 0.25d0 (+ pn lam)))))))))
  1931.  
  1932. (defun check-drj (x y z p)
  1933.   (let* ((w (/ (* x y) z))
  1934.      (a (* p p (+ x y z w)))
  1935.      (b (* p (+ p x) (+ p y)))
  1936.      (drj-1 (drj x (+ x z) (+ x w) (+ x p)))
  1937.      (drj-2 (drj y (+ y z) (+ y w) (+ y p)))
  1938.      (drj-3 (drj a b b a))
  1939.      (drj-4 (drj 0d0 z w p)))
  1940.     ;; Both values should be equal.
  1941.     (values (+ drj-1 drj-2 (* (- a b) drj-3) (/ 3 (sqrt a)))
  1942.         drj-4)))
  1943.  
  1944. ;;; Other Jacobian elliptic functions
  1945.  
  1946. ;; jacobi_ns(u,m) = 1/jacobi_sn(u,m)
  1947. (defmfun  $jacobi_ns (u m)
  1948.   (simplify (list '(%jacobi_ns) (resimplify u) (resimplify m))))
  1949.  
  1950. (defprop %jacobi_ns simp-%jacobi_ns operators)
  1951.  
  1952. (defprop %jacobi_ns
  1953.     ((u m)
  1954.      ;; diff wrt u
  1955.      ((mtimes) -1 ((%jacobi_cn) u m) ((%jacobi_dn) u m)
  1956.       ((mexpt) ((%jacobi_sn) u m) -2))
  1957.      ;; diff wrt m
  1958.      ((mtimes) -1 ((mexpt) ((%jacobi_sn) u m) -2)
  1959.       ((mplus)
  1960.        ((mtimes) ((rat) 1 2)
  1961.     ((mexpt) ((mplus) 1 ((mtimes) -1 m)) -1)
  1962.     ((mexpt) ((%jacobi_cn) u m) 2)
  1963.     ((%jacobi_sn) u m))
  1964.        ((mtimes) ((rat) 1 2) ((mexpt) m -1)
  1965.     ((%jacobi_cn) u m) ((%jacobi_dn) u m)
  1966.     ((mplus) u
  1967.      ((mtimes) -1
  1968.       ((mexpt) ((mplus) 1 ((mtimes) -1 m)) -1)
  1969.       (($elliptic_e) ((%asin) ((%jacobi_sn) u m))
  1970.        m)))))))
  1971.   grad)
  1972.  
  1973. (defmfun simp-%jacobi_ns (form y z)
  1974.   (declare (ignore y))
  1975.   (twoargcheck form)
  1976.   (let ((u (simpcheck (cadr form) z))
  1977.     (m (simpcheck (caddr form) z))
  1978.     coef)
  1979.     (cond ((or (and (floatp u) (floatp m))
  1980.            (and $numer (numberp u) (numberp m)))
  1981.        ;; Numerically evaluate sn
  1982.        (/ (sn (float u 1d0) (float m 1d0))))
  1983.       ((and $numer (complex-number-p u)
  1984.         (complex-number-p m))
  1985.        (let ((u-r ($realpart u))
  1986.          (u-i ($imagpart u))
  1987.          (m-r ($realpart m))
  1988.          (m-i ($imagpart m)))
  1989.          (complexify (/ (sn (complex u-r u-i) (complex m-r m-i))))))
  1990.       ((zerop1 m)
  1991.        ;; A&S 16.6.10
  1992.        `(($csc) ,u))
  1993.       ((onep1 m)
  1994.        ;; A&S 16.6.10
  1995.        `(($coth) ,u))
  1996.       ((zerop1 u)
  1997.        (dbz-err1 'jacobi_ns))
  1998.       ((and $trigsign (mminusp* u))
  1999.        ;; ns is odd
  2000.        (neg (cons-exp '%jacobi_ns (neg u) m)))
  2001.       ;; A&S 16.20 (Jacobi's Imaginary transformation)
  2002.       ((and $%iargs (multiplep u '$%i))
  2003.        ;; ns(i*u) = 1/sn(i*u) = -i/sc(u,m1) = -i*cs(u,m1)
  2004.        (neg (mul '$%i
  2005.              (cons-exp '%jacobi_cs (coeff u '$%i 1) (add 1 (neg m))))))
  2006.       ((setq coef (kc-arg2 u m))
  2007.        ;; A&S 16.8.10
  2008.        ;;
  2009.        ;; ns(m*K+u) = 1/sn(m*K+u)
  2010.        ;;
  2011.        (destructuring-bind (lin const)
  2012.            coef
  2013.          (cond ((integerp lin)
  2014.             (ecase (mod lin 4)
  2015.               (0
  2016.                ;; ns(4*m*K+u) = ns(u)
  2017.                ;; ns(0) = infinity
  2018.                (if (zerop1 const)
  2019.                (dbz-err1 'jacobi_ns)
  2020.                `((%jacobi_ns simp) ,const ,m)))
  2021.               (1
  2022.                ;; ns(4*m*K + K + u) = ns(K+u) = dc(u)
  2023.                ;; ns(K) = 1
  2024.                (if (zerop1 const)
  2025.                1
  2026.                `((%jacobi_dc simp) ,const ,m)))
  2027.               (2
  2028.                ;; ns(4*m*K+2*K + u) = ns(2*K+u) = -ns(u)
  2029.                ;; ns(2*K) = infinity
  2030.                (if (zerop1 const)
  2031.                (dbz-err1 'jacobi_ns)
  2032.                (neg `((%jacobi_ns simp) ,const ,m))))
  2033.               (3
  2034.                ;; ns(4*m*K+3*K+u) = ns(2*K + K + u) = -ns(K+u) = -dc(u)
  2035.                ;; ns(3*K) = -1
  2036.                (if (zerop1 const)
  2037.                -1
  2038.                (neg `((%jacobi_dc simp) ,const ,m))))))
  2039.            ((and (alike1 lin 1//2)
  2040.              (zerop1 const))
  2041.             `((mexpt) ((%jacobi_sn) ,u ,m) -1))
  2042.            (t
  2043.             (eqtest (list '(%jacobi_ns) u m) form)))))      
  2044.       (t
  2045.        ;; Nothing to do
  2046.        (eqtest (list '(%jacobi_ns) u m) form)))))
  2047.  
  2048. ;; jacobi_nc(u,m) = 1/jacobi_cn(u,m)
  2049.  
  2050. (defmfun  $jacobi_nc (u m)
  2051.   (simplify (list '(%jacobi_nc) (resimplify u) (resimplify m))))
  2052.  
  2053. (defprop %jacobi_nc simp-%jacobi_nc operators)
  2054.  
  2055. (defprop %jacobi_nc
  2056.     ((u m)
  2057.      ;; wrt u
  2058.      ((mtimes) ((mexpt) ((%jacobi_cn) u m) -2)
  2059.       ((%jacobi_dn) u m) ((%jacobi_sn) u m))
  2060.      ;; wrt m
  2061.      ((mtimes) -1 ((mexpt) ((%jacobi_cn) u m) -2)
  2062.       ((mplus)
  2063.        ((mtimes) ((rat) -1 2)
  2064.     ((mexpt) ((mplus) 1 ((mtimes) -1 m)) -1)
  2065.     ((%jacobi_cn) u m) ((mexpt) ((%jacobi_sn) u m) 2))
  2066.        ((mtimes) ((rat) -1 2) ((mexpt) m -1)
  2067.     ((%jacobi_dn) u m) ((%jacobi_sn) u m)
  2068.     ((mplus) u
  2069.      ((mtimes) -1 ((mexpt) ((mplus) 1 ((mtimes) -1 m)) -1)
  2070.       (($elliptic_e) ((%asin) ((%jacobi_sn) u m)) m)))))))
  2071.   grad)
  2072.  
  2073. (defmfun simp-%jacobi_nc (form y z)
  2074.   (declare (ignore y))
  2075.   (twoargcheck form)
  2076.   (let ((u (simpcheck (cadr form) z))
  2077.     (m (simpcheck (caddr form) z))
  2078.     coef)
  2079.     (cond ((or (and (floatp u) (floatp m))
  2080.            (and $numer (numberp u) (numberp m)))
  2081.        (/ (cn (float u 1d0) (float m 1d0))))
  2082.       ((and $numer (complex-number-p u)
  2083.         (complex-number-p m))
  2084.        (let ((u-r ($realpart u))
  2085.          (u-i ($imagpart u))
  2086.          (m-r ($realpart m))
  2087.          (m-i ($imagpart m)))
  2088.          (complexify (/ (cn (complex u-r u-i) (complex m-r m-i))))))
  2089.       ((zerop1 u)
  2090.        1)
  2091.       ((zerop1 m)
  2092.        ;; A&S 16.6.8
  2093.        `(($sec) ,u))
  2094.       ((onep1 m)
  2095.        ;; A&S 16.6.8
  2096.        `((%cosh) ,u))
  2097.       ((and $trigsign (mminusp* u))
  2098.        ;; nc is even
  2099.        (cons-exp '%jacobi_nc (neg u) m))
  2100.       ;; A&S 16.20 (Jacobi's Imaginary transformation)
  2101.       ((and $%iargs (multiplep u '$%i))
  2102.        ;; nc(i*u) = 1/cn(i*u) = 1/nc(u,1-m) = cn(u,1-m)
  2103.        (cons-exp '%jacobi_cn (coeff u '$%i 1) (add 1 (neg m))))
  2104.       ((setq coef (kc-arg2 u m))
  2105.        ;; A&S 16.8.8
  2106.        ;;
  2107.        ;; nc(u) = 1/cn(u)
  2108.        ;;
  2109.        (destructuring-bind (lin const)
  2110.            coef
  2111.          (cond ((integerp lin)
  2112.             (ecase (mod lin 4)
  2113.               (0
  2114.                ;; nc(4*m*K+u) = nc(u)
  2115.                ;; nc(0) = 1
  2116.                (if (zerop1 const)
  2117.                1
  2118.                `((%jacobi_nc simp) ,const ,m)))
  2119.               (1
  2120.                ;; nc(4*m*K+K+u) = nc(K+u) = -ds(u)/sqrt(1-m)
  2121.                ;; nc(K) = infinity
  2122.                (if (zerop1 const)
  2123.                (dbz-err1 'jacobi_nc)
  2124.                (neg `((mtimes simp)
  2125.                   ((mexpt simp)
  2126.                    ((mplus simp) 1 ((mtimes simp) -1 ,m))
  2127.                    ((rat simp) -1 2))
  2128.                   ((%jacobi_ds simp) ,const ,m)))))
  2129.               (2
  2130.                ;; nc(4*m*K+2*K+u) = nc(2*K+u) = -nc(u)
  2131.                ;; nc(2*K) = -1
  2132.                (if (zerop1 const)
  2133.                -1
  2134.                (neg `((%jacobi_nc) ,const ,m))))
  2135.               (3
  2136.                ;; nc(4*m*K+3*K+u) = nc(3*K+u) = nc(2*K+K+u) =
  2137.                ;; -nc(K+u) = ds(u)/sqrt(1-m)
  2138.                ;;
  2139.                ;; nc(3*K) = infinity
  2140.                (if (zerop1 const)
  2141.                (dbz-err1 'jacobi_nc)
  2142.                `((mtimes simp)
  2143.                  ((mexpt simp)
  2144.                   ((mplus simp) 1 ((mtimes simp) -1 ,m))
  2145.                   ((rat simp) -1 2))
  2146.                  ((%jacobi_ds simp) ,const ,m))))))
  2147.            ((and (alike1 1//2 lin)
  2148.              (zerop1 const))
  2149.             `((mexpt) ((%jacobi_cn) ,u ,m) -1))
  2150.            (t
  2151.             (eqtest (list '(%jacobi_cn) u m) form)))))
  2152.       (t
  2153.        ;; Nothing to do
  2154.        (eqtest (list '(%jacobi_nc) u m) form)))))
  2155.  
  2156. ;; jacobi_nd(u,m) = 1/jacobi_dn(u,m)
  2157. (defmfun  $jacobi_nd (u m)
  2158.   (simplify (list '(%jacobi_nd) (resimplify u) (resimplify m))))
  2159.  
  2160. (defprop %jacobi_nd simp-%jacobi_nd operators)
  2161.  
  2162. (defprop %jacobi_nd
  2163.     ((u m)
  2164.      ;; wrt u
  2165.      ((mtimes) m ((%jacobi_cn) u m)
  2166.       ((mexpt) ((%jacobi_dn) u m) -2) ((%jacobi_sn) u m))
  2167.      ;; wrt m
  2168.      ((mtimes) -1 ((mexpt) ((%jacobi_dn) u m) -2)
  2169.       ((mplus)
  2170.        ((mtimes) ((rat) -1 2)
  2171.     ((mexpt) ((mplus) 1 ((mtimes) -1 m)) -1)
  2172.     ((%jacobi_dn) u m)
  2173.     ((mexpt) ((%jacobi_sn) u m) 2))
  2174.        ((mtimes) ((rat) -1 2) ((%jacobi_cn) u m)
  2175.     ((%jacobi_sn) u m)
  2176.     ((mplus) u
  2177.      ((mtimes) -1
  2178.       ((mexpt) ((mplus) 1 ((mtimes) -1 m)) -1)
  2179.       (($elliptic_e) ((%asin) ((%jacobi_sn) u m))
  2180.        m)))))))
  2181.   grad)
  2182.  
  2183. (defmfun simp-%jacobi_nd (form y z)
  2184.   (declare (ignore y))
  2185.   (twoargcheck form)
  2186.   (let ((u (simpcheck (cadr form) z))
  2187.     (m (simpcheck (caddr form) z))
  2188.     coef)
  2189.     (cond ((or (and (floatp u) (floatp m))
  2190.            (and $numer (numberp u) (numberp m)))
  2191.        (/ (dn (float u 1d0) (float m 1d0))))
  2192.       ((and $numer (complex-number-p u)
  2193.         (complex-number-p m))
  2194.        (let ((u-r ($realpart u))
  2195.          (u-i ($imagpart u))
  2196.          (m-r ($realpart m))
  2197.          (m-i ($imagpart m)))
  2198.          (complexify (/ (dn (complex u-r u-i) (complex m-r m-i))))))
  2199.       ((zerop1 u)
  2200.        1)
  2201.       ((zerop1 m)
  2202.        ;; A&S 16.6.6
  2203.        1)
  2204.       ((onep1 m)
  2205.        ;; A&S 16.6.6
  2206.        `((%cosh) ,u))
  2207.       ((and $trigsign (mminusp* u))
  2208.        ;; nd is even
  2209.        (cons-exp '%jacobi_nd (neg u) m))
  2210.       ;; A&S 16.20 (Jacobi's Imaginary transformation)
  2211.       ((and $%iargs (multiplep u '$%i))
  2212.        ;; nd(i*u) = 1/dn(i*u) = 1/dc(u,1-m) = cd(u,1-m)
  2213.        (cons-exp '%jacobi_cd (coeff u '$%i 1) (add 1 (neg m))))
  2214.       ((setq coef (kc-arg2 u m))
  2215.        ;; A&S 16.8.6
  2216.        ;;
  2217.        (destructuring-bind (lin const)
  2218.            coef
  2219.          (cond ((integerp lin)
  2220.             ;; nd has period 2K
  2221.             (ecase (mod lin 2)
  2222.               (0
  2223.                ;; nd(2*m*K+u) = nd(u)
  2224.                ;; nd(0) = 1
  2225.                (if (zerop1 const)
  2226.                1
  2227.                `((%jacobi_nd) ,const ,m)))
  2228.               (1
  2229.                ;; nd(2*m*K+K+u) = nd(K+u) = dn(u)/sqrt(1-m)
  2230.                ;; nd(K) = 1/sqrt(1-m)
  2231.                (if (zerop1 const)
  2232.                `((mexpt simp)
  2233.                  ((mplus simp) 1 ((mtimes simp) -1 ,m))
  2234.                  ((rat simp) -1 2))
  2235.                `((mtimes simp)
  2236.                  ((%jacobi_nd simp) ,const ,m)
  2237.                  ((mexpt simp)
  2238.                   ((mplus simp) 1 ((mtimes simp) -1 ,m))
  2239.                   ((rat simp) -1 2)))))))
  2240.            (t
  2241.             (eqtest (list '(%jacobi_nd) u m) form)))))
  2242.       (t
  2243.        ;; Nothing to do
  2244.        (eqtest (list '(%jacobi_nd) u m) form)))))
  2245.  
  2246. ;; jacobi_sc(u,m) = jacobi_sn/jacobi_cn
  2247. (defmfun  $jacobi_sc (u m)
  2248.   (simplify (list '(%jacobi_sc) (resimplify u) (resimplify m))))
  2249.  
  2250. (defprop %jacobi_sc simp-%jacobi_sc operators)
  2251.  
  2252. (defprop %jacobi_sc
  2253.     ((u m)
  2254.      ;; wrt u
  2255.      ((mtimes) ((mexpt) ((%jacobi_cn) u m) -2)
  2256.       ((%jacobi_dn) u m))
  2257.      ;; wrt m
  2258.      ((mplus)
  2259.       ((mtimes) ((mexpt) ((%jacobi_cn) u m) -1)
  2260.        ((mplus)
  2261.     ((mtimes) ((rat) 1 2)
  2262.      ((mexpt) ((mplus) 1 ((mtimes) -1 m)) -1)
  2263.      ((mexpt) ((%jacobi_cn) u m) 2)
  2264.      ((%jacobi_sn) u m))
  2265.     ((mtimes) ((rat) 1 2) ((mexpt) m -1)
  2266.      ((%jacobi_cn) u m) ((%jacobi_dn) u m)
  2267.      ((mplus) u
  2268.       ((mtimes) -1
  2269.        ((mexpt) ((mplus) 1 ((mtimes) -1 m)) -1)
  2270.        (($elliptic_e) ((%asin) ((%jacobi_sn) u m))
  2271.         m))))))
  2272.       ((mtimes) -1 ((mexpt) ((%jacobi_cn) u m) -2)
  2273.        ((%jacobi_sn) u m)
  2274.        ((mplus)
  2275.     ((mtimes) ((rat) -1 2)
  2276.      ((mexpt) ((mplus) 1 ((mtimes) -1 m)) -1)
  2277.      ((%jacobi_cn) u m)
  2278.      ((mexpt) ((%jacobi_sn) u m) 2))
  2279.     ((mtimes) ((rat) -1 2) ((mexpt) m -1)
  2280.      ((%jacobi_dn) u m) ((%jacobi_sn) u m)
  2281.      ((mplus) u
  2282.       ((mtimes) -1
  2283.        ((mexpt) ((mplus) 1 ((mtimes) -1 m)) -1)
  2284.        (($elliptic_e) ((%asin) ((%jacobi_sn) u m))
  2285.         m))))))))
  2286.   grad)
  2287.  
  2288. (defmfun simp-%jacobi_sc (form y z)
  2289.   (declare (ignore y))
  2290.   (twoargcheck form)
  2291.   (let ((u (simpcheck (cadr form) z))
  2292.     (m (simpcheck (caddr form) z))
  2293.     coef)
  2294.     (cond ((or (and (floatp u) (floatp m))
  2295.            (and $numer (numberp u) (numberp m)))
  2296.        (let ((fu (float u 1d0))
  2297.          (fm (float m 1d0)))
  2298.          (/ (sn fu fm) (cn fu fm))))
  2299.       ((and $numer (complex-number-p u)
  2300.         (complex-number-p m))
  2301.        (let ((u-r ($realpart u))
  2302.          (u-i ($imagpart u))
  2303.          (m-r ($realpart m))
  2304.          (m-i ($imagpart m)))
  2305.          (complexify (/ (sn (complex u-r u-i) (complex m-r m-i))
  2306.                 (cn (complex u-r u-i) (complex m-r m-i))))))
  2307.       ((zerop1 u)
  2308.        0)
  2309.       ((zerop1 m)
  2310.        ;; A&S 16.6.9
  2311.        `((%tan) ,u))
  2312.       ((onep1 m)
  2313.        ;; A&S 16.6.9
  2314.        `((%sinh) ,u))
  2315.       ((and $trigsign (mminusp* u))
  2316.        ;; sc is odd
  2317.        (neg (cons-exp '%jacobi_sc (neg u) m)))
  2318.       ;; A&S 16.20 (Jacobi's Imaginary transformation)
  2319.       ((and $%iargs (multiplep u '$%i))
  2320.        ;; sc(i*u) = sn(i*u)/cn(i*u) = i*sc(u,m1)/nc(u,m1) = i*sn(u,m1)
  2321.        (mul '$%i
  2322.         (cons-exp '%jacobi_sn (coeff u '$%i 1) (add 1 (neg m)))))
  2323.       ((setq coef (kc-arg2 u m))
  2324.        ;; A&S 16.8.9
  2325.        ;; sc(2*m*K+u) = sc(u)
  2326.        (destructuring-bind (lin const)
  2327.            coef
  2328.          (cond ((integerp lin)
  2329.             (ecase (mod lin 2)
  2330.               (0
  2331.                ;; sc(2*m*K+ u) = sc(u)
  2332.                ;; sc(0) = 0
  2333.                (if (zerop1 const)
  2334.                1
  2335.                `((%jacobi_sc simp) ,const ,m)))
  2336.               (1
  2337.                ;; sc(2*m*K + K + u) = sc(K+u)= - cs(u)/sqrt(1-m)
  2338.                ;; sc(K) = infinity
  2339.                (if (zerop1 const)
  2340.                (dbz-err1 'jacobi_sc)
  2341.                `((mtimes simp) -1
  2342.                  ((mexpt simp)
  2343.                   ((mplus simp) 1 ((mtimes simp) -1 ,m))
  2344.                   ((rat simp) -1 2))
  2345.                  ((%jacobi_cs simp) ,const ,m))))))
  2346.            ((and (alike1 lin 1//2)
  2347.              (zerop1 const))
  2348.             ;; (1-m)^(1/4)
  2349.             `((mexpt simp)
  2350.               ((mplus simp) 1 ((mtimes simp) -1 ,m))
  2351.               ((rat simp) 1 4)))
  2352.            (t
  2353.             (eqtest (list '(%jacobi_sc) u m) form)))))
  2354.       (t
  2355.        ;; Nothing to do
  2356.        (eqtest (list '(%jacobi_sc) u m) form)))))
  2357.  
  2358. ;; jacobi_sd(u,m) = jacobi_sn/jacobi_dn
  2359. (defmfun  $jacobi_sd (u m)
  2360.   (simplify (list '(%jacobi_sd) (resimplify u) (resimplify m))))
  2361.  
  2362. (defprop %jacobi_sd simp-%jacobi_sd operators)
  2363.  
  2364. (defprop %jacobi_sd
  2365.     ((u m)
  2366.      ;; wrt u
  2367.      ((mtimes) ((%jacobi_cn) u m)
  2368.       ((mexpt) ((%jacobi_dn) u m) -2))
  2369.      ;; wrt m
  2370.      ((mplus)
  2371.       ((mtimes) ((mexpt) ((%jacobi_dn) u m) -1)
  2372.        ((mplus)
  2373.     ((mtimes) ((rat) 1 2)
  2374.      ((mexpt) ((mplus) 1 ((mtimes) -1 m)) -1)
  2375.      ((mexpt) ((%jacobi_cn) u m) 2)
  2376.      ((%jacobi_sn) u m))
  2377.     ((mtimes) ((rat) 1 2) ((mexpt) m -1)
  2378.      ((%jacobi_cn) u m) ((%jacobi_dn) u m)
  2379.      ((mplus) u
  2380.       ((mtimes) -1
  2381.        ((mexpt) ((mplus) 1 ((mtimes) -1 m)) -1)
  2382.        (($elliptic_e) ((%asin) ((%jacobi_sn) u m))
  2383.         m))))))
  2384.       ((mtimes) -1 ((mexpt) ((%jacobi_dn) u m) -2)
  2385.        ((%jacobi_sn) u m)
  2386.        ((mplus)
  2387.     ((mtimes) ((rat) -1 2)
  2388.      ((mexpt) ((mplus) 1 ((mtimes) -1 m)) -1)
  2389.      ((%jacobi_dn) u m)
  2390.      ((mexpt) ((%jacobi_sn) u m) 2))
  2391.     ((mtimes) ((rat) -1 2) ((%jacobi_cn) u m)
  2392.      ((%jacobi_sn) u m)
  2393.      ((mplus) u
  2394.       ((mtimes) -1
  2395.        ((mexpt) ((mplus) 1 ((mtimes) -1 m)) -1)
  2396.        (($elliptic_e) ((%asin) ((%jacobi_sn) u m))
  2397.         m))))))))
  2398.   grad)
  2399.  
  2400. (defmfun simp-%jacobi_sd (form y z)
  2401.   (declare (ignore y))
  2402.   (twoargcheck form)
  2403.   (let ((u (simpcheck (cadr form) z))
  2404.     (m (simpcheck (caddr form) z))
  2405.     coef)
  2406.     (cond ((or (and (floatp u) (floatp m))
  2407.            (and $numer (numberp u) (numberp m)))
  2408.        (let ((fu (float u 1d0))
  2409.          (fm (float m 1d0)))
  2410.          (/ (sn fu fm) (dn fu fm))))
  2411.       ((and $numer (complex-number-p u)
  2412.         (complex-number-p m))
  2413.        (let ((u-r ($realpart u))
  2414.          (u-i ($imagpart u))
  2415.          (m-r ($realpart m))
  2416.          (m-i ($imagpart m)))
  2417.          (complexify (/ (sn (complex u-r u-i) (complex m-r m-i))
  2418.                 (dn (complex u-r u-i) (complex m-r m-i))))))
  2419.       ((zerop1 u)
  2420.        0)
  2421.       ((zerop1 m)
  2422.        ;; A&S 16.6.5
  2423.        `((%sin) ,u))
  2424.       ((onep1 m)
  2425.        ;; A&S 16.6.5
  2426.        `((%sinh) ,u))
  2427.       ((and $trigsign (mminusp* u))
  2428.        ;; sd is odd
  2429.        (neg (cons-exp '%jacobi_sd (neg u) m)))
  2430.       ;; A&S 16.20 (Jacobi's Imaginary transformation)
  2431.       ((and $%iargs (multiplep u '$%i))
  2432.        ;; sd(i*u) = sn(i*u)/dn(i*u) = i*sc(u,m1)/dc(u,m1) = i*sd(u,m1)
  2433.        (mul '$%i
  2434.         (cons-exp '%jacobi_sd (coeff u '$%i 1) (add 1 (neg m)))))
  2435.       ((setq coef (kc-arg2 u m))
  2436.        ;; A&S 16.8.5
  2437.        ;; sd(4*m*K+u) = sd(u)
  2438.        (destructuring-bind (lin const)
  2439.            coef
  2440.          (cond ((integerp lin)
  2441.             (ecase (mod lin 4)
  2442.               (0
  2443.                ;; sd(4*m*K+u) = sd(u)
  2444.                ;; sd(0) = 0
  2445.                (if (zerop1 const)
  2446.                0
  2447.                `((%jacobi_sd simp) ,const ,m)))
  2448.               (1
  2449.                ;; sd(4*m*K+K+u) = sd(K+u) = cn(u)/sqrt(1-m)
  2450.                ;; sd(K) = 1/sqrt(m1)
  2451.                (if (zerop1 const)
  2452.                `((mexpt) ((mplus) 1 ((mtimes) -1 ,m))
  2453.                  ((rat) -1 2))
  2454.                `((mtimes simp)
  2455.                  ((mexpt simp)
  2456.                   ((mplus simp) 1 ((mtimes simp) -1 ,m))
  2457.                   ((rat simp) -1 2))
  2458.                  ((%jacobi_cn simp) ,const ,m))))
  2459.               (2
  2460.                ;; sd(4*m*K+2*K+u) = sd(2*K+u) = -sd(u)
  2461.                ;; sd(2*K) = 0
  2462.                (if (zerop1 const)
  2463.                0
  2464.                (neg `((%jacobi_sd) ,const ,m))))
  2465.               (3
  2466.                ;; sd(4*m*K+3*K+u) = sd(3*K+u) = sd(2*K+K+u) =
  2467.                ;; -sd(K+u) = -cn(u)/sqrt(1-m)
  2468.                ;; sd(3*K) = -1/sqrt(m1)
  2469.                (if (zerop1 const)
  2470.                (neg `((mexpt)
  2471.                   ((mplus simp) 1 ((mtimes simp) -1 ,m))
  2472.                   ((rat) -1 2)))
  2473.                (neg `((mtimes simp)
  2474.                   ((mexpt simp)
  2475.                    ((mplus simp) 1 ((mtimes simp) -1 ,m))
  2476.                    ((rat simp) -1 2))
  2477.                   ((%jacobi_cn simp) ,const ,m)))))))
  2478.            ((and (alike1 lin 1//2)
  2479.              (zerop1 const))
  2480.             ;; jacobi_sn/jacobi_dn
  2481.             `((mtimes)
  2482.               ((%jacobi_sn) ((mtimes) ((rat) 1 2)
  2483.                      ((%elliptic_kc) ,m))
  2484.                ,m)
  2485.               ((mexpt)
  2486.                ((%jacobi_dn) ((mtimes) ((rat) 1 2)
  2487.                       ((%elliptic_kc) ,m))
  2488.             ,m)
  2489.                -1)))
  2490.            (t
  2491.             (eqtest (list '(%jacobi_sd) u m) form)))))
  2492.       (t
  2493.        ;; Nothing to do
  2494.        (eqtest (list '(%jacobi_sd) u m) form)))))
  2495.  
  2496. ;; jacobi_cs(u,m) = jacobi_cn/jacobi_sn
  2497. (defmfun  $jacobi_cs (u m)
  2498.   (simplify (list '(%jacobi_cs) (resimplify u) (resimplify m))))
  2499.  
  2500. (defprop %jacobi_cs simp-%jacobi_cs operators)
  2501.  
  2502. (defprop %jacobi_cs
  2503.     ((u m)
  2504.      ;; wrt u
  2505.      ((mtimes) -1 ((%jacobi_dn) u m)
  2506.       ((mexpt) ((%jacobi_sn) u m) -2))
  2507.      ;; wrt m
  2508.      ((mplus)
  2509.  ((mtimes) -1 ((%jacobi_cn) u m)
  2510.   ((mexpt) ((%jacobi_sn) u m) -2)
  2511.   ((mplus)
  2512.    ((mtimes) ((rat) 1 2)
  2513.     ((mexpt) ((mplus) 1 ((mtimes) -1 m)) -1)
  2514.     ((mexpt) ((%jacobi_cn) u m) 2)
  2515.     ((%jacobi_sn) u m))
  2516.    ((mtimes) ((rat) 1 2) ((mexpt) m -1)
  2517.     ((%jacobi_cn) u m) ((%jacobi_dn) u m)
  2518.     ((mplus) u
  2519.      ((mtimes) -1
  2520.       ((mexpt) ((mplus) 1 ((mtimes) -1 m)) -1)
  2521.       (($elliptic_e) ((%asin) ((%jacobi_sn) u m))
  2522.        m))))))
  2523.  ((mtimes) ((mexpt) ((%jacobi_sn) u m) -1)
  2524.   ((mplus)
  2525.    ((mtimes) ((rat) -1 2)
  2526.     ((mexpt) ((mplus) 1 ((mtimes) -1 m)) -1)
  2527.     ((%jacobi_cn) u m)
  2528.     ((mexpt) ((%jacobi_sn) u m) 2))
  2529.    ((mtimes) ((rat) -1 2) ((mexpt) m -1)
  2530.     ((%jacobi_dn) u m) ((%jacobi_sn) u m)
  2531.     ((mplus) u
  2532.      ((mtimes) -1
  2533.       ((mexpt) ((mplus) 1 ((mtimes) -1 m)) -1)
  2534.       (($elliptic_e) ((%asin) ((%jacobi_sn) u m))
  2535.        m))))))))
  2536.   grad)
  2537.  
  2538. (defmfun simp-%jacobi_cs (form y z)
  2539.   (declare (ignore y))
  2540.   (twoargcheck form)
  2541.   (let ((u (simpcheck (cadr form) z))
  2542.     (m (simpcheck (caddr form) z))
  2543.     coef)
  2544.     (cond ((or (and (floatp u) (floatp m))
  2545.            (and $numer (numberp u) (numberp m)))
  2546.        (let ((fu (float u 1d0))
  2547.          (fm (float m 1d0)))
  2548.          (/ (cn fu fm) (sn fu fm))))
  2549.       ((and $numer (complex-number-p u)
  2550.         (complex-number-p m))
  2551.        (let ((u-r ($realpart u))
  2552.          (u-i ($imagpart u))
  2553.          (m-r ($realpart m))
  2554.          (m-i ($imagpart m)))
  2555.          (complexify (/ (cn (complex u-r u-i) (complex m-r m-i))
  2556.                 (sn (complex u-r u-i) (complex m-r m-i))))))
  2557.       ((zerop1 m)
  2558.        ;; A&S 16.6.12
  2559.        `(($cot) ,u))
  2560.       ((onep1 m)
  2561.        ;; A&S 16.6.12
  2562.        `(($csch) ,u))
  2563.       ((zerop1 u)
  2564.        (dbz-err1 'jacobi_cs))
  2565.       ((and $trigsign (mminusp* u))
  2566.        ;; cs is odd
  2567.        (neg (cons-exp '%jacobi_cs (neg u) m)))
  2568.       ;; A&S 16.20 (Jacobi's Imaginary transformation)
  2569.       ((and $%iargs (multiplep u '$%i))
  2570.        ;; cs(i*u) = cn(i*u)/sn(i*u) = -i*nc(u,m1)/sc(u,m1) = -i*ns(u,m1)
  2571.        (neg (mul '$%i
  2572.              (cons-exp '%jacobi_ns (coeff u '$%i 1) (add 1 (neg m))))))
  2573.       ((setq coef (kc-arg2 u m))
  2574.        ;; A&S 16.8.12
  2575.        ;; 
  2576.        ;; cs(2*m*K + u) = cs(u)
  2577.        (destructuring-bind (lin const)
  2578.            coef
  2579.          (cond ((integerp lin)
  2580.             (ecase (mod lin 2)
  2581.               (0
  2582.                ;; cs(2*m*K + u) = cs(u)
  2583.                ;; cs(0) = infinity
  2584.                (if (zerop1 const)
  2585.                (dbz-err1 'jacobi_cs)
  2586.                `((%jacobi_cs simp) ,const ,m)))
  2587.               (1
  2588.                ;; cs(K+u) = -sqrt(1-m)*sc(u)
  2589.                ;; cs(K) = 0
  2590.                (if (zerop1 const)
  2591.                0
  2592.                `((mtimes simp) -1
  2593.                  ((mexpt simp)
  2594.                   ((mplus simp) 1 ((mtimes simp) -1 ,m))
  2595.                   ((rat simp) 1 2))
  2596.                  ((%jacobi_sc simp) ,const ,m))))))
  2597.            ((and (alike1 lin 1//2)
  2598.              (zerop1 const))
  2599.             ;; 1/jacobi_sc
  2600.             `((mexpt)
  2601.               ((%jacobi_sc) ((mtimes) ((rat) 1 2)
  2602.                      ((%elliptic_kc) ,m)) ,m)
  2603.               -1))
  2604.            (t
  2605.             (eqtest (list '(%jacobi_cs simp) u m) form)))))
  2606.       (t
  2607.        ;; Nothing to do
  2608.        (eqtest (list '(%jacobi_cs simp) u m) form)))))
  2609.  
  2610. ;; jacobi_cd(u,m) = jacobi_cn/jacobi_dn
  2611. (defmfun  $jacobi_cd (u m)
  2612.   (simplify (list '(%jacobi_cd) (resimplify u) (resimplify m))))
  2613.  
  2614. (defprop %jacobi_cd simp-%jacobi_cd operators)
  2615.  
  2616. (defprop %jacobi_cd
  2617.     ((u m)
  2618.      ;; wrt u
  2619.      ((mtimes) ((mplus) -1 m)
  2620.       ((mexpt) ((%jacobi_dn) u m) -2)
  2621.       ((%jacobi_sn) u m))
  2622.      ;; wrt m
  2623.      ((mplus)
  2624.       ((mtimes) -1 ((%jacobi_cn) u m)
  2625.        ((mexpt) ((%jacobi_dn) u m) -2)
  2626.        ((mplus)
  2627.     ((mtimes) ((rat) -1 2)
  2628.      ((mexpt) ((mplus) 1 ((mtimes) -1 m)) -1)
  2629.      ((%jacobi_dn) u m)
  2630.      ((mexpt) ((%jacobi_sn) u m) 2))
  2631.     ((mtimes) ((rat) -1 2) ((%jacobi_cn) u m)
  2632.      ((%jacobi_sn) u m)
  2633.      ((mplus) u
  2634.       ((mtimes) -1
  2635.        ((mexpt) ((mplus) 1 ((mtimes) -1 m)) -1)
  2636.        (($elliptic_e) ((%asin) ((%jacobi_sn) u m))
  2637.         m))))))
  2638.       ((mtimes) ((mexpt) ((%jacobi_dn) u m) -1)
  2639.        ((mplus)
  2640.     ((mtimes) ((rat) -1 2)
  2641.      ((mexpt) ((mplus) 1 ((mtimes) -1 m)) -1)
  2642.      ((%jacobi_cn) u m)
  2643.      ((mexpt) ((%jacobi_sn) u m) 2))
  2644.     ((mtimes) ((rat) -1 2) ((mexpt) m -1)
  2645.      ((%jacobi_dn) u m) ((%jacobi_sn) u m)
  2646.      ((mplus) u
  2647.       ((mtimes) -1
  2648.        ((mexpt) ((mplus) 1 ((mtimes) -1 m)) -1)
  2649.        (($elliptic_e) ((%asin) ((%jacobi_sn) u m))
  2650.         m))))))))
  2651.   grad)
  2652.  
  2653. (defmfun simp-%jacobi_cd (form y z)
  2654.   (declare (ignore y))
  2655.   (twoargcheck form)
  2656.   (let ((u (simpcheck (cadr form) z))
  2657.     (m (simpcheck (caddr form) z))
  2658.     coef)
  2659.     (cond ((or (and (floatp u) (floatp m))
  2660.            (and $numer (numberp u) (numberp m)))
  2661.        (let ((fu (float u 1d0))
  2662.          (fm (float m 1d0)))
  2663.          (/ (cn fu fm) (dn fu fm))))
  2664.       ((and $numer (complex-number-p u)
  2665.         (complex-number-p m))
  2666.        (let ((u-r ($realpart u))
  2667.          (u-i ($imagpart u))
  2668.          (m-r ($realpart m))
  2669.          (m-i ($imagpart m)))
  2670.          (complexify (/ (cn (complex u-r u-i) (complex m-r m-i))
  2671.                 (dn (complex u-r u-i) (complex m-r m-i))))))
  2672.       ((zerop1 u)
  2673.        1)
  2674.       ((zerop1 m)
  2675.        ;; A&S 16.6.4
  2676.        `((%cos) ,u))
  2677.       ((onep1 m)
  2678.        ;; A&S 16.6.4
  2679.        1)
  2680.       ((and $trigsign (mminusp* u))
  2681.        ;; cd is even
  2682.        (cons-exp '%jacobi_cd (neg u) m))
  2683.       ;; A&S 16.20 (Jacobi's Imaginary transformation)
  2684.       ((and $%iargs (multiplep u '$%i))
  2685.        ;; cd(i*u) = cn(i*u)/dn(i*u) = nc(u,m1)/dc(u,m1) = nd(u,m1)
  2686.        (cons-exp '%jacobi_nd (coeff u '$%i 1) (add 1 (neg m))))
  2687.       ((setf coef (kc-arg2 u m))
  2688.        ;; A&S 16.8.4
  2689.        ;;
  2690.        (destructuring-bind (lin const)
  2691.            coef
  2692.          (cond ((integerp lin)
  2693.             (ecase (mod lin 4)
  2694.               (0
  2695.                ;; cd(4*m*K + u) = cd(u)
  2696.                ;; cd(0) = 1
  2697.                (if (zerop1 const)
  2698.                1
  2699.                `((%jacobi_cd) ,const ,m)))
  2700.               (1
  2701.                ;; cd(4*m*K + K + u) = cd(K+u) = -sn(u)
  2702.                ;; cd(K) = 0
  2703.                (if (zerop1 const)
  2704.                0
  2705.                (neg `((%jacobi_sn) ,const ,m))))
  2706.               (2
  2707.                ;; cd(4*m*K + 2*K + u) = cd(2*K+u) = -cd(u)
  2708.                ;; cd(2*K) = -1
  2709.                (if (zerop1 const)
  2710.                -1
  2711.                (neg `((%jacobi_cd) ,const ,m))))
  2712.               (3
  2713.                ;; cd(4*m*K + 3*K + u) = cd(2*K + K + u) =
  2714.                ;; -cd(K+u) = sn(u)
  2715.                ;; cd(3*K) = 0
  2716.                (if (zerop1 const)
  2717.                0
  2718.                `((%jacobi_sn) ,const ,m)))))
  2719.            ((and (alike1 lin 1//2)
  2720.              (zerop1 const))
  2721.             ;; jacobi_cn/jacobi_dn
  2722.             `((mtimes)
  2723.               ((%jacobi_cn) ((mtimes) ((rat) 1 2)
  2724.                      ((%elliptic_kc) ,m))
  2725.                ,m)
  2726.               ((mexpt)
  2727.                ((%jacobi_dn) ((mtimes) ((rat) 1 2)
  2728.                       ((%elliptic_kc) ,m))
  2729.             ,m)
  2730.                -1)))
  2731.            (t
  2732.             ;; Nothing to do
  2733.             (eqtest (list '(%jacobi_cd) u m) form)))))
  2734.       (t
  2735.        ;; Nothing to do
  2736.        (eqtest (list '(%jacobi_cd) u m) form)))))
  2737.  
  2738. ;; jacobi_ds(u,m) = jacobi_dn/jacobi_sn
  2739. (defmfun  $jacobi_ds (u m)
  2740.   (simplify (list '(%jacobi_ds) (resimplify u) (resimplify m))))
  2741.  
  2742. (defprop %jacobi_ds simp-%jacobi_ds operators)
  2743.  
  2744. (defprop %jacobi_ds
  2745.     ((u m)
  2746.      ;; wrt u
  2747.      ((mtimes) -1 ((%jacobi_cn) u m)
  2748.       ((mexpt) ((%jacobi_sn) u m) -2))
  2749.      ;; wrt m
  2750.      ((mplus)
  2751.       ((mtimes) -1 ((%jacobi_dn) u m)
  2752.        ((mexpt) ((%jacobi_sn) u m) -2)
  2753.        ((mplus)
  2754.     ((mtimes) ((rat) 1 2)
  2755.      ((mexpt) ((mplus) 1 ((mtimes) -1 m)) -1)
  2756.      ((mexpt) ((%jacobi_cn) u m) 2)
  2757.      ((%jacobi_sn) u m))
  2758.     ((mtimes) ((rat) 1 2) ((mexpt) m -1)
  2759.      ((%jacobi_cn) u m) ((%jacobi_dn) u m)
  2760.      ((mplus) u
  2761.       ((mtimes) -1
  2762.        ((mexpt) ((mplus) 1 ((mtimes) -1 m)) -1)
  2763.        (($elliptic_e) ((%asin) ((%jacobi_sn) u m))
  2764.         m))))))
  2765.       ((mtimes) ((mexpt) ((%jacobi_sn) u m) -1)
  2766.        ((mplus)
  2767.     ((mtimes) ((rat) -1 2)
  2768.      ((mexpt) ((mplus) 1 ((mtimes) -1 m)) -1)
  2769.      ((%jacobi_dn) u m)
  2770.      ((mexpt) ((%jacobi_sn) u m) 2))
  2771.     ((mtimes) ((rat) -1 2) ((%jacobi_cn) u m)
  2772.      ((%jacobi_sn) u m)
  2773.      ((mplus) u
  2774.       ((mtimes) -1
  2775.        ((mexpt) ((mplus) 1 ((mtimes) -1 m)) -1)
  2776.        (($elliptic_e) ((%asin) ((%jacobi_sn) u m))
  2777.         m))))))))
  2778.   grad)
  2779.  
  2780. (defmfun simp-%jacobi_ds (form y z)
  2781.   (declare (ignore y))
  2782.   (twoargcheck form)
  2783.   (let ((u (simpcheck (cadr form) z))
  2784.     (m (simpcheck (caddr form) z))
  2785.     coef)
  2786.     (cond ((or (and (floatp u) (floatp m))
  2787.            (and $numer (numberp u) (numberp m)))
  2788.        (let ((fu (float u 1d0))
  2789.          (fm (float m 1d0)))
  2790.          (/ (dn fu fm) (sn fu fm))))
  2791.       ((and $numer (complex-number-p u)
  2792.         (complex-number-p m))
  2793.        (let ((u-r ($realpart u))
  2794.          (u-i ($imagpart u))
  2795.          (m-r ($realpart m))
  2796.          (m-i ($imagpart m)))
  2797.          (complexify (/ (dn (complex u-r u-i) (complex m-r m-i))
  2798.                 (sn (complex u-r u-i) (complex m-r m-i))))))
  2799.       ((zerop1 m)
  2800.        ;; A&S 16.6.11
  2801.        `(($csc) ,u))
  2802.       ((onep1 m)
  2803.        ;; A&S 16.6.11
  2804.        `(($csch) ,u))
  2805.       ((zerop1 u)
  2806.        (dbz-err1 'jacobi_ds))
  2807.       ((and $trigsign (mminusp* u))
  2808.        (neg (cons-exp '%jacobi_ds (neg u) m)))
  2809.       ;; A&S 16.20 (Jacobi's Imaginary transformation)
  2810.       ((and $%iargs (multiplep u '$%i))
  2811.        ;; ds(i*u) = dn(i*u)/sn(i*u) = -i*dc(u,m1)/sc(u,m1) = -i*ds(u,m1)
  2812.        (neg (mul '$%i
  2813.              (cons-exp '%jacobi_ds (coeff u '$%i 1) (add 1 (neg m))))))
  2814.       ((setf coef (kc-arg2 u m))
  2815.        ;; A&S 16.8.11
  2816.        (destructuring-bind (lin const)
  2817.            coef
  2818.          (cond ((integerp lin)
  2819.             (ecase (mod lin 4)
  2820.               (0
  2821.                ;; ds(4*m*K + u) = ds(u)
  2822.                ;; ds(0) = infinity
  2823.                (if (zerop1 const)
  2824.                (dbz-err1 'jacobi_ds)
  2825.                `((%jacobi_ds) ,const ,m)))
  2826.               (1
  2827.                ;; ds(4*m*K + K + u) = ds(K+u) = sqrt(1-m)*nc(u)
  2828.                ;; ds(K) = sqrt(1-m)
  2829.                (if (zerop1 const)
  2830.                `((mexpt simp)
  2831.                  ((mplus simp) 1 ((mtimes simp) -1 ,m))
  2832.                  ((rat simp) 1 2))
  2833.                `((mtimes simp)
  2834.                  ((mexpt simp)
  2835.                   ((mplus simp) 1 ((mtimes simp) -1 ,m))
  2836.                   ((rat simp) 1 2))
  2837.                  ((%jacobi_nc simp) ,const ,m))))
  2838.               (2
  2839.                ;; ds(4*m*K + 2*K + u) = ds(2*K+u) = -ds(u)
  2840.                ;; ds(0) = pole
  2841.                (if (zerop1 const)
  2842.                (dbz-err1 'jacobi_ds)
  2843.                (neg `((%jacobi_ds) ,const ,m))))
  2844.               (3
  2845.                ;; ds(4*m*K + 3*K + u) = ds(2*K + K + u) =
  2846.                ;; -ds(K+u) = -sqrt(1-m)*nc(u)
  2847.                ;; ds(3*K) = -sqrt(1-m)
  2848.                (if (zerop1 const)
  2849.                (neg `((mexpt simp)
  2850.                   ((mplus simp) 1 ((mtimes simp) -1 ,m))
  2851.                   ((rat simp) 1 2)))
  2852.                (neg `((mtimes simp)
  2853.                   ((mexpt simp)
  2854.                    ((mplus simp) 1 ((mtimes simp) -1 ,m))
  2855.                    ((rat simp) 1 2))
  2856.                   ((%jacobi_nc simp) ,const ,m)))))))
  2857.            ((and (alike1 lin 1//2)
  2858.              (zerop1 const))
  2859.             ;; jacobi_dn/jacobi_sn
  2860.             `((mtimes)
  2861.               ((%jacobi_dn) ((mtimes) ((rat) 1 2)
  2862.                      ((%elliptic_kc) ,m))
  2863.                ,m)
  2864.               ((mexpt)
  2865.                ((%jacobi_sn) ((mtimes) ((rat) 1 2)
  2866.                       ((%elliptic_kc) ,m))
  2867.             ,m)
  2868.                -1)))
  2869.            (t
  2870.             ;; Nothing to do
  2871.             (eqtest (list '(%jacobi_ds) u m) form)))))
  2872.       (t
  2873.        ;; Nothing to do
  2874.        (eqtest (list '(%jacobi_ds) u m) form)))))
  2875.  
  2876. ;; jacobi_dc(u,m) = jacobi_dn/jacobi_cn
  2877. (defmfun  $jacobi_dc (u m)
  2878.   (simplify (list '(%jacobi_dc) (resimplify u) (resimplify m))))
  2879.  
  2880. (defprop %jacobi_dc simp-%jacobi_dc operators)
  2881.  
  2882. (defprop %jacobi_dc
  2883.     ((u m)
  2884.      ;; wrt u
  2885.      ((mtimes) ((mplus) 1 ((mtimes) -1 m))
  2886.       ((mexpt) ((%jacobi_cn) u m) -2)
  2887.       ((%jacobi_sn) u m))
  2888.      ;; wrt m
  2889.      ((mplus)
  2890.       ((mtimes) ((mexpt) ((%jacobi_cn) u m) -1)
  2891.        ((mplus)
  2892.     ((mtimes) ((rat) -1 2)
  2893.      ((mexpt) ((mplus) 1 ((mtimes) -1 m)) -1)
  2894.      ((%jacobi_dn) u m)
  2895.      ((mexpt) ((%jacobi_sn) u m) 2))
  2896.     ((mtimes) ((rat) -1 2) ((%jacobi_cn) u m)
  2897.      ((%jacobi_sn) u m)
  2898.      ((mplus) u
  2899.       ((mtimes) -1
  2900.        ((mexpt) ((mplus) 1 ((mtimes) -1 m)) -1)
  2901.        (($elliptic_e) ((%asin) ((%jacobi_sn) u m))
  2902.         m))))))
  2903.       ((mtimes) -1 ((mexpt) ((%jacobi_cn) u m) -2)
  2904.        ((%jacobi_dn) u m)
  2905.        ((mplus)
  2906.     ((mtimes) ((rat) -1 2)
  2907.      ((mexpt) ((mplus) 1 ((mtimes) -1 m)) -1)
  2908.      ((%jacobi_cn) u m)
  2909.      ((mexpt) ((%jacobi_sn) u m) 2))
  2910.     ((mtimes) ((rat) -1 2) ((mexpt) m -1)
  2911.      ((%jacobi_dn) u m) ((%jacobi_sn) u m)
  2912.      ((mplus) u
  2913.       ((mtimes) -1
  2914.        ((mexpt) ((mplus) 1 ((mtimes) -1 m)) -1)
  2915.        (($elliptic_e) ((%asin) ((%jacobi_sn) u m))
  2916.         m))))))))
  2917.   grad)
  2918.  
  2919. (defmfun simp-%jacobi_dc (form y z)
  2920.   (declare (ignore y))
  2921.   (twoargcheck form)
  2922.   (let ((u (simpcheck (cadr form) z))
  2923.     (m (simpcheck (caddr form) z))
  2924.     coef)
  2925.     (cond ((or (and (floatp u) (floatp m))
  2926.            (and $numer (numberp u) (numberp m)))
  2927.        (let ((fu (float u 1d0))
  2928.          (fm (float m 1d0)))
  2929.          (/ (dn fu fm) (cn fu fm))))
  2930.       ((zerop1 u)
  2931.        1)
  2932.       ((zerop1 m)
  2933.        ;; A&S 16.6.7
  2934.        `(($sec) ,u))
  2935.       ((onep1 m)
  2936.        ;; A&S 16.6.7
  2937.        1)
  2938.       ((and $trigsign (mminusp* u))
  2939.        (cons-exp '%jacobi_dc (neg u) m))
  2940.       ;; A&S 16.20 (Jacobi's Imaginary transformation)
  2941.       ((and $%iargs (multiplep u '$%i))
  2942.        ;; dc(i*u) = dn(i*u)/cn(i*u) = dc(u,m1)/nc(u,m1) = dn(u,m1)
  2943.        (cons-exp '%jacobi_dn (coeff u '$%i 1) (add 1 (neg m))))
  2944.       ((setf coef (kc-arg2 u m))
  2945.        ;; See A&S 16.8.7
  2946.        (destructuring-bind (lin const)
  2947.            coef
  2948.          (cond ((integerp lin)
  2949.             (ecase (mod lin 4)
  2950.               (0
  2951.                ;; dc(4*m*K + u) = dc(u)
  2952.                ;; dc(0) = 1
  2953.                (if (zerop1 const)
  2954.                1
  2955.                `((%jacobi_dc) ,const ,m)))
  2956.               (1
  2957.                ;; dc(4*m*K + K + u) = dc(K+u) = -ns(u)
  2958.                ;; dc(K) = pole
  2959.                (if (zerop1 const)
  2960.                (dbz-err1 'jacobi_dc)
  2961.                (neg `((%jacobi_ns simp) ,const ,m))))
  2962.               (2
  2963.                ;; dc(4*m*K + 2*K + u) = dc(2*K+u) = -dc(u)
  2964.                ;; dc(2K) = -1
  2965.                (if (zerop1 const)
  2966.                -1
  2967.                (neg `((%jacobi_dc) ,const ,m))))
  2968.               (3
  2969.                ;; dc(4*m*K + 3*K + u) = dc(2*K + K + u) =
  2970.                ;; -dc(K+u) = ns(u)
  2971.                ;; dc(3*K) = ns(0) = inf
  2972.                (if (zerop1 const)
  2973.                (dbz-err1 'jacobi_dc)
  2974.                `((%jacobi_dc simp) ,const ,m)))))
  2975.            ((and (alike1 lin 1//2)
  2976.              (zerop1 const))
  2977.             ;; jacobi_dn/jacobi_cn
  2978.             `((mtimes)
  2979.               ((%jacobi_dn) ((mtimes) ((rat) 1 2)
  2980.                      ((%elliptic_kc) ,m))
  2981.                ,m)
  2982.               ((mexpt)
  2983.                ((%jacobi_cn) ((mtimes) ((rat) 1 2)
  2984.                       ((%elliptic_kc) ,m))
  2985.             ,m)
  2986.                -1)))
  2987.            (t
  2988.             ;; Nothing to do
  2989.             (eqtest (list '(%jacobi_dc) u m) form)))))
  2990.       (t
  2991.        ;; Nothing to do
  2992.        (eqtest (list '(%jacobi_dc) u m) form)))))
  2993.  
  2994. ;;; Other inverse Jacobian functions
  2995.  
  2996. ;; inverse_jacobi_ns(x)
  2997. ;;
  2998. ;; Let u = inverse_jacobi_ns(x).  Then jacobi_ns(u) = x or
  2999. ;; 1/jacobi_sn(u) = x or
  3000. ;;
  3001. ;; jacobi_sn(u) = 1/x
  3002. ;;
  3003. ;; so u = inverse_jacobi_sn(1/x)
  3004.  
  3005. (defmfun $inverse_jacobi_ns (u m)
  3006.   (simplify (list '(%inverse_jacobi_ns) (resimplify u) (resimplify m))))
  3007.  
  3008. (defprop %inverse_jacobi_ns
  3009.     ((x m)
  3010.      ;; -1/sqrt(1-x^2)/sqrt(x^2+m-1)
  3011.      ((mtimes) -1
  3012.       ((mexpt) ((mplus) -1 ((mexpt) x 2)) ((rat) -1 2))
  3013.       ((mexpt)
  3014.        ((mplus) ((mtimes simp ratsimp) -1 m) ((mexpt) x 2))
  3015.        ((rat) -1 2)))
  3016.      ;; wrt m
  3017.      ((%derivative) ((%inverse_jacobi_ns) x m) m 1))
  3018.   grad)
  3019.  
  3020. (defprop %inverse_jacobi_ns simp-%inverse_jacobi_ns operators)
  3021.  
  3022. (defmfun simp-%inverse_jacobi_ns (form y z)
  3023.   (declare (ignore y))
  3024.   (twoargcheck form)
  3025.   (let ((u (simpcheck (cadr form) z))
  3026.     (m (simpcheck (caddr form) z)))
  3027.     (cond ((or (and (floatp u) (floatp m))
  3028.            (and $numer (numberp u) (numberp m)))
  3029.        ;; Numerically evaluate asn
  3030.        ;;
  3031.        ;; ans(x,m) = asn(1/x,m) = F(asin(1/x),m)
  3032.        (elliptic-f (lisp:asin (/ u)) m))
  3033.       ((zerop1 m)
  3034.        ;; ans(x,0) = F(asin(1/x),0) = asin(1/x)
  3035.        `((%elliptic_f) ((%asin) ((mexpt) ,u -1)) 0))
  3036.       ((onep1 m)
  3037.        ;; ans(x,1) = F(asin(1/x),1) = log(tan(pi/2+asin(1/x)/2))
  3038.        `((%elliptic_f) ((%asin) ((mexpt) ,u -1)) 1))
  3039.       ((onep1 u)
  3040.        `((%elliptic_kc) ,m))
  3041.       ((alike1 u -1)
  3042.        (neg `((%elliptic_kc) ,m)))
  3043.       (t
  3044.        ;; Nothing to do
  3045.        (eqtest (list '(%inverse_jacobi_ns) u m) form)))))
  3046.  
  3047. ;; inverse_jacobi_nc(x)
  3048. ;;
  3049. ;; Let u = inverse_jacobi_nc(x).  Then jacobi_nc(u) = x or
  3050. ;; 1/jacobi_cn(u) = x or
  3051. ;;
  3052. ;; jacobi_cn(u) = 1/x
  3053. ;;
  3054. ;; so u = inverse_jacobi_cn(1/x)
  3055.  
  3056. (defmfun $inverse_jacobi_nc (u m)
  3057.   (simplify (list '(%inverse_jacobi_nc) (resimplify u) (resimplify m))))
  3058.  
  3059. (defprop %inverse_jacobi_nc
  3060.     ((x m)
  3061.      ;; -1/sqrt(1-x^2)/sqrt(x^2+m-1)
  3062.      ((mtimes)
  3063.       ((mexpt) ((mplus) -1 ((mexpt) x 2)) ((rat) -1 2))
  3064.       ((mexpt)
  3065.        ((mplus) ((mtimes simp ratsimp) -1 m) ((mexpt) x 2))
  3066.        ((rat) -1 2)))
  3067.      ;; wrt m
  3068.      ((%derivative) ((%inverse_jacobi_nc) x m) m 1))
  3069.   grad)
  3070.  
  3071. (defprop %inverse_jacobi_nc simp-%inverse_jacobi_nc operators)
  3072.  
  3073. (defmfun simp-%inverse_jacobi_nc (form y z)
  3074.   (declare (ignore y))
  3075.   (twoargcheck form)
  3076.   (let ((u (simpcheck (cadr form) z))
  3077.     (m (simpcheck (caddr form) z)))
  3078.     (cond ((or (and (floatp u) (floatp m))
  3079.            (and $numer (numberp u) (numberp m)))
  3080.        ;;
  3081.        ($inverse_jacobi_cn (/ u) m))
  3082.       ((onep1 u)
  3083.        0)
  3084.       ((alike1 u -1)
  3085.        `((mtimes) 2 ((%elliptic_kc) ,m)))
  3086.       (t
  3087.        ;; Nothing to do
  3088.        (eqtest (list '(%inverse_jacobi_nc) u m) form)))))
  3089.  
  3090. ;; inverse_jacobi_nd(x)
  3091. ;;
  3092. ;; Let u = inverse_jacobi_nd(x).  Then jacobi_nd(u) = x or
  3093. ;; 1/jacobi_dn(u) = x or
  3094. ;;
  3095. ;; jacobi_dn(u) = 1/x
  3096. ;;
  3097. ;; so u = inverse_jacobi_dn(1/x)
  3098.  
  3099. (defmfun $inverse_jacobi_nd (u m)
  3100.   (simplify (list '(%inverse_jacobi_nd) (resimplify u) (resimplify m))))
  3101.  
  3102. (defprop %inverse_jacobi_nd
  3103.     ((x m)
  3104.      ;; -1/sqrt(1-x^2)/sqrt(x^2+m-1)
  3105.      ((mtimes) -1
  3106.       ((mexpt) ((mplus) -1 ((mexpt simp ratsimp) x 2))
  3107.        ((rat) -1 2))
  3108.       ((mexpt)
  3109.        ((mplus) 1
  3110.     ((mtimes) ((mplus) -1 m) ((mexpt simp ratsimp) x 2)))
  3111.        ((rat) -1 2)))
  3112.      ;; wrt m
  3113.      ((%derivative) ((%inverse_jacobi_nd) x m) m 1))
  3114.   grad)
  3115.  
  3116. (defprop %inverse_jacobi_nd simp-%inverse_jacobi_nd operators)
  3117.  
  3118. (defmfun simp-%inverse_jacobi_nd (form y z)
  3119.   (declare (ignore y))
  3120.   (twoargcheck form)
  3121.   (let ((u (simpcheck (cadr form) z))
  3122.     (m (simpcheck (caddr form) z)))
  3123.     (cond ((or (and (floatp u) (floatp m))
  3124.            (and $numer (numberp u) (numberp m)))
  3125.        ($inverse_jacobi_dn (/ u) m))
  3126.       ((onep1 u)
  3127.        `((%elliptic_kc) ,m))
  3128.       (t
  3129.        ;; Nothing to do
  3130.        (eqtest (list '(%inverse_jacobi_nd) u m) form)))))
  3131.  
  3132. ;; inverse_jacobi_sc(x)
  3133. ;;
  3134. ;; Let u = inverse_jacobi_sc(x).  Then jacobi_sc(u) = x or
  3135. ;; x = jacobi_sn(u)/jacobi_cn(u)
  3136. ;;
  3137. ;; x^2 = sn^2/cn^2
  3138. ;;     = sn^2/(1-sn^2)
  3139. ;;
  3140. ;; so
  3141. ;;
  3142. ;; sn^2 = x^2/(1+x^2)
  3143. ;;
  3144. ;; sn(u) = x/sqrt(1+x^2)
  3145. ;;
  3146. ;; u = inverse_sn(x/sqrt(1+x^2))
  3147. ;;
  3148.  
  3149. (defmfun $inverse_jacobi_sc (u m)
  3150.   (simplify (list '(%inverse_jacobi_sc) (resimplify u) (resimplify m))))
  3151.  
  3152. (defprop %inverse_jacobi_sc
  3153.     ((x m)
  3154.      ((mtimes)
  3155.       ((mexpt) ((mplus) 1 ((mexpt) x 2))
  3156.        ((rat) -1 2))
  3157.       ((mexpt)
  3158.        ((mplus) 1
  3159.     ((mtimes) -1 ((mplus) -1 m) ((mexpt) x 2)))
  3160.        ((rat) -1 2)))
  3161.      ;; wrt m
  3162.      ((%derivative) ((%inverse_jacobi_sc) x m) m 1))
  3163.   grad)
  3164.  
  3165. (defprop %inverse_jacobi_sc simp-%inverse_jacobi_sc operators)
  3166.  
  3167. (defmfun simp-%inverse_jacobi_sc (form y z)
  3168.   (declare (ignore y))
  3169.   (twoargcheck form)
  3170.   (let ((u (simpcheck (cadr form) z))
  3171.     (m (simpcheck (caddr form) z)))
  3172.     (cond ((or (and (floatp u) (floatp m))
  3173.            (and $numer (numberp u) (numberp m)))
  3174.        ($inverse_jacobi_sn (/ u (sqrt (+ 1 (* u u)))) m))
  3175.       (t
  3176.        ;; Nothing to do
  3177.        (eqtest (list '(%inverse_jacobi_sc) u m) form)))))
  3178.  
  3179. ;; inverse_jacobi_sd(x)
  3180. ;;
  3181. ;; Let u = inverse_jacobi_sd(x).  Then jacobi_sd(u) = x or
  3182. ;; x = jacobi_sn(u)/jacobi_dn(u)
  3183. ;;
  3184. ;; x^2 = sn^2/dn^2
  3185. ;;     = sn^2/(1-m*sn^2)
  3186. ;;
  3187. ;; so
  3188. ;;
  3189. ;; sn^2 = x^2/(1+m*x^2)
  3190. ;;
  3191. ;; sn(u) = x/sqrt(1+m*x^2)
  3192. ;;
  3193. ;; u = inverse_sn(x/sqrt(1+m*x^2))
  3194. ;;
  3195.  
  3196. (defmfun $inverse_jacobi_sd (u m)
  3197.   (simplify (list '(%inverse_jacobi_sd) (resimplify u) (resimplify m))))
  3198.  
  3199. (defprop %inverse_jacobi_sd
  3200.     ((x m)
  3201.      ((mtimes)
  3202.       ((mexpt)
  3203.        ((mplus) 1 ((mtimes) ((mplus) -1 m) ((mexpt) x 2)))
  3204.        ((rat) -1 2))
  3205.       ((mexpt) ((mplus) 1 ((mtimes) m ((mexpt) x 2)))
  3206.        ((rat) -1 2)))
  3207.      ;; wrt m
  3208.      ((%derivative) ((%inverse_jacobi_sd) x m) m 1))
  3209.   grad)
  3210.  
  3211. (defprop %inverse_jacobi_sd simp-%inverse_jacobi_sd operators)
  3212.  
  3213. (defmfun simp-%inverse_jacobi_sd (form y z)
  3214.   (declare (ignore y))
  3215.   (twoargcheck form)
  3216.   (let ((u (simpcheck (cadr form) z))
  3217.     (m (simpcheck (caddr form) z)))
  3218.     (cond ((or (and (floatp u) (floatp m))
  3219.            (and $numer (numberp u) (numberp m)))
  3220.        ($inverse_jacobi_sn (/ u (sqrt (+ 1 (* m u u)))) m))
  3221.       ((zerop1 u)
  3222.        0)
  3223.       (t
  3224.        ;; Nothing to do
  3225.        (eqtest (list '(%inverse_jacobi_sd) u m) form)))))
  3226.  
  3227. ;; inverse_jacobi_cs(x)
  3228. ;;
  3229. ;; Let u = inverse_jacobi_cs(x).  Then jacobi_cs(u) = x or
  3230. ;; 1/x = 1/jacobi_cs(u) = jacobi_sc(u)
  3231. ;;
  3232. ;; u = inverse_sc(1/x)
  3233. ;;
  3234.  
  3235. (defmfun $inverse_jacobi_cs (u m)
  3236.   (simplify (list '(%inverse_jacobi_cs) (resimplify u) (resimplify m))))
  3237.  
  3238. (defprop %inverse_jacobi_cs
  3239.     ((x m)
  3240.      ((mtimes) -1
  3241.       ((mexpt) ((mplus) 1 ((mexpt simp ratsimp) x 2))
  3242.        ((rat) -1 2))
  3243.       ((mexpt) ((mplus) 1
  3244.              ((mtimes simp ratsimp) -1 m)
  3245.              ((mexpt simp ratsimp) x 2))
  3246.        ((rat) -1 2)))
  3247.      ;; wrt m
  3248.      ((%derivative) ((%inverse_jacobi_cs) x m) m 1))
  3249.   grad)
  3250.  
  3251. (defprop %inverse_jacobi_cs simp-%inverse_jacobi_cs operators)
  3252.  
  3253. (defmfun simp-%inverse_jacobi_cs (form y z)
  3254.   (declare (ignore y))
  3255.   (twoargcheck form)
  3256.   (let ((u (simpcheck (cadr form) z))
  3257.     (m (simpcheck (caddr form) z)))
  3258.     (cond ((or (and (floatp u) (floatp m))
  3259.            (and $numer (numberp u) (numberp m)))
  3260.        ($inverse_jacobi_sc (/ u) m))
  3261.       ((zerop1 u)
  3262.        `((%elliptic_kc) ,m))
  3263.       (t
  3264.        ;; Nothing to do
  3265.        (eqtest (list '(%inverse_jacobi_cs) u m) form)))))
  3266.  
  3267. ;; inverse_jacobi_cd(x)
  3268. ;;
  3269. ;; Let u = inverse_jacobi_cd(x).  Then jacobi_cd(u) = x or
  3270. ;; x = jacobi_cn(u)/jacobi_dn(u)
  3271. ;;
  3272. ;; x^2 = cn^2/dn^2
  3273. ;;     = (1-sn^2)/(1-m*sn^2)
  3274. ;;
  3275. ;; or
  3276. ;;
  3277. ;; sn^2 = (1-x^2)/(1-m*x^2)
  3278. ;;
  3279. ;; sn(u) = sqrt(1-x^2)/sqrt(1-m*x^2)
  3280. ;;
  3281. ;; u = inverse_sn(sqrt(1-x^2)/sqrt(1-m*x^2))
  3282. ;;
  3283.  
  3284. (defmfun $inverse_jacobi_cd (u m)
  3285.   (simplify (list '(%inverse_jacobi_cd) (resimplify u) (resimplify m))))
  3286.  
  3287. (defprop %inverse_jacobi_cd
  3288.     ((x m)
  3289.      ((mtimes)
  3290.       ((mexpt)
  3291.        ((mplus) 1 ((mtimes) -1 ((mexpt) x 2)))
  3292.        ((rat) -1 2))
  3293.       ((mexpt)
  3294.        ((mplus) 1 ((mtimes) -1 m ((mexpt) x 2)))
  3295.        ((rat) -1 2)))
  3296.      ;; wrt m
  3297.      ((%derivative) ((%inverse_jacobi_cd) x m) m 1))
  3298.   grad)
  3299.  
  3300. (defprop %inverse_jacobi_cd simp-%inverse_jacobi_cd operators)
  3301.  
  3302. (defmfun simp-%inverse_jacobi_cd (form y z)
  3303.   (declare (ignore y))
  3304.   (twoargcheck form)
  3305.   (let ((u (simpcheck (cadr form) z))
  3306.     (m (simpcheck (caddr form) z)))
  3307.     (cond ((or (and (floatp u) (floatp m))
  3308.            (and $numer (numberp u) (numberp m)))
  3309.        ($inverse_jacobi_sn (/ (sqrt (* (- 1 u) (+ 1 u)))
  3310.                   (sqrt (- 1 (* m u u)))) m))
  3311.       ((onep1 u)
  3312.        0)
  3313.       ((zerop1 u)
  3314.        `((%elliptic_kc) ,m))
  3315.       (t
  3316.        ;; Nothing to do
  3317.        (eqtest (list '(%inverse_jacobi_cd) u m) form)))))
  3318.  
  3319. ;; inverse_jacobi_ds(x)
  3320. ;;
  3321. ;; Let u = inverse_jacobi_ds(x).  Then jacobi_ds(u) = x or
  3322. ;; 1/x = 1/jacobi_ds(u) = jacobi_sd(u)
  3323. ;;
  3324. ;; u = inverse_sd(1/x)
  3325. ;;
  3326.  
  3327. (defmfun $inverse_jacobi_ds (u m)
  3328.   (simplify (list '(%inverse_jacobi_ds) (resimplify u) (resimplify m))))
  3329.  
  3330. (defprop %inverse_jacobi_ds
  3331.     ((x m)
  3332.      ((mtimes) -1
  3333.       ((mexpt)
  3334.        ((mplus) -1 m ((mexpt simp ratsimp) x 2))
  3335.        ((rat) -1 2))
  3336.       ((mexpt)
  3337.        ((mplus) m ((mexpt simp ratsimp) x 2))
  3338.        ((rat) -1 2)))
  3339.       ;; wrt m
  3340.      ((%derivative) ((%inverse_jacobi_ds) x m) m 1))
  3341.   grad)
  3342.  
  3343. (defprop %inverse_jacobi_ds simp-%inverse_jacobi_ds operators)
  3344.  
  3345. (defmfun simp-%inverse_jacobi_ds (form y z)
  3346.   (declare (ignore y))
  3347.   (twoargcheck form)
  3348.   (let ((u (simpcheck (cadr form) z))
  3349.     (m (simpcheck (caddr form) z)))
  3350.     (cond ((or (and (floatp u) (floatp m))
  3351.            (and $numer (numberp u) (numberp m)))
  3352.        ($inverse_jacobi_sd (/ u) m))
  3353.       ((and $trigsign (mminusp* u))
  3354.        (neg (cons-exp '%inverse_jacobi_ds (neg u) m)))
  3355.       (t
  3356.        ;; Nothing to do
  3357.        (eqtest (list '(%inverse_jacobi_ds) u m) form)))))
  3358.  
  3359.  
  3360. ;; inverse_jacobi_dc(x)
  3361. ;;
  3362. ;; Let u = inverse_jacobi_dc(x).  Then jacobi_dc(u) = x or
  3363. ;; 1/x = 1/jacobi_dc(u) = jacobi_cd(u)
  3364. ;;
  3365. ;; u = inverse_cd(1/x)
  3366. ;;
  3367.  
  3368. (defmfun $inverse_jacobi_dc (u m)
  3369.   (simplify (list '(%inverse_jacobi_dc) (resimplify u) (resimplify m))))
  3370.  
  3371. (defprop %inverse_jacobi_dc
  3372.     ((x m)
  3373.      ((mtimes) -1
  3374.       ((mexpt)
  3375.        ((mplus) -1 ((mexpt simp ratsimp) x 2))
  3376.        ((rat) -1 2))
  3377.       ((mexpt)
  3378.        ((mplus)
  3379.     ((mtimes simp ratsimp) -1 m)
  3380.     ((mexpt simp ratsimp) x 2))
  3381.        ((rat) -1 2)))
  3382.      ;; wrt m
  3383.      ((%derivative) ((%inverse_jacobi_dc) x m) m 1))
  3384.   grad)
  3385.  
  3386. (defprop %inverse_jacobi_dc simp-%inverse_jacobi_dc operators)
  3387.  
  3388. (defmfun simp-%inverse_jacobi_dc (form y z)
  3389.   (declare (ignore y))
  3390.   (twoargcheck form)
  3391.   (let ((u (simpcheck (cadr form) z))
  3392.     (m (simpcheck (caddr form) z)))
  3393.     (cond ((or (and (floatp u) (floatp m))
  3394.            (and $numer (numberp u) (numberp m)))
  3395.        ($inverse_jacobi_cd (/ u) m))
  3396.       (t
  3397.        ;; Nothing to do
  3398.        (eqtest (list '(%inverse_jacobi_dc) u m) form)))))
  3399.  
  3400. ;; Convert an inverse Jacobian function into the equivalent elliptic
  3401. ;; integral F.
  3402. ;;
  3403. ;; See A&S 17.4.41-17.4.52.
  3404. (defun make-elliptic-f (e)
  3405.   (cond ((atom e)
  3406.      e)
  3407.     ((member (caar e) '(%inverse_jacobi_sc %inverse_jacobi_cs
  3408.                 %inverse_jacobi_nd %inverse_jacobi_dn
  3409.                 %inverse_jacobi_sn %inverse_jacobi_cd
  3410.                 %inverse_jacobi_dc %inverse_jacobi_ns
  3411.                 %inverse_jacobi_nc %inverse_jacobi_ds
  3412.                 %inverse_jacobi_sd %inverse_jacobi_cn))
  3413.      ;; We have some inverse Jacobi function.  Convert it to the F form.
  3414.      (destructuring-bind ((fn &rest ops) u m)
  3415.          e
  3416.        (declare (ignore ops))
  3417.        (ecase fn
  3418.          (%inverse_jacobi_sc
  3419.           ;; A&S 17.4.41
  3420.           `(($elliptic_f) ((%atan) ,u) ,m))
  3421.          (%inverse_jacobi_cs
  3422.           ;; A&S 17.4.42
  3423.           `(($elliptic_f) ((%atan) ((mexpt) ,u -1)) ,m))
  3424.          (%inverse_jacobi_nd
  3425.           ;; A&S 17.4.43
  3426.           `(($elliptic_f)
  3427.         ((%asin) ((mtimes)
  3428.               ((mexpt) ,m ((rat) -1 2))
  3429.               ((mexpt) ,u -1)
  3430.               ((mexpt) ((mplus) -1 ((mexpt) ,u 2))
  3431.                ((rat) 1 2))))
  3432.         ,m))
  3433.          (%inverse_jacobi_dn
  3434.           ;; A&S 17.4.44
  3435.           `(($elliptic_f)
  3436.         ((%asin)
  3437.          ((mtimes) ((mexpt) ,m ((rat) -1 2))
  3438.           ((mexpt) ((mplus) 1 ((mtimes) -1 ((mexpt) ,u 2)))
  3439.            ((rat) 1 2))))
  3440.         ,m))
  3441.          (%inverse_jacobi_sn
  3442.           ;; A&S 17.4.45
  3443.           `(($elliptic_f) ((%asin) ,u) ,m))
  3444.          (%inverse_jacobi_cd
  3445.           ;; A&S 17.4.46
  3446.           `(($elliptic_f)
  3447.         ((%asin)
  3448.          ((mexpt) ((mtimes) ((mplus) 1
  3449.                      ((mtimes) -1 ((mexpt) ,u 2)))
  3450.                ((mexpt) ((mplus) 1
  3451.                      ((mtimes) -1 ,m ((mexpt) ,u 2)))
  3452.                 -1))
  3453.           ((rat) 1 2)))
  3454.         ,m))
  3455.          (%inverse_jacobi_dc
  3456.           ;; A&S 17.4.47
  3457.           `(($elliptic_f)
  3458.         ((%asin)
  3459.          ((mexpt)
  3460.           ((mtimes) ((mplus) -1 ((mexpt) ,u 2))
  3461.            ((mexpt) ((mplus) ((mtimes) -1 ,m) ((mexpt) ,u 2)) -1))
  3462.           ((rat) 1 2)))
  3463.         ,m))
  3464.          (%inverse_jacobi_ns
  3465.           ;; A&S 17.4.48
  3466.           `(($elliptic_f) ((asin) ((mexpt) ,u -1)) ,m))
  3467.          (%inverse_jacobi_nc
  3468.           ;; A&S 17.4.49
  3469.           `(($elliptic_f) ((acos) ((mexpt) ,u -1)) ,m))
  3470.          (%inverse_jacobi_ds
  3471.           ;; A&S 17.4.50
  3472.           `(($elliptic_f)
  3473.         ((%asin) ((mexpt) ((mplus) ,m ((mexpt) ,u 2))
  3474.               ((rat) -1 2)))
  3475.         ,m))
  3476.          (%inverse_jacobi_sd
  3477.           ;; A&S 17.4.51
  3478.           `(($elliptic_f)
  3479.         ((%asin)
  3480.          ((mtimes) ,u
  3481.           ((mexpt) ((mplus) 1 ((mtimes) ,m ((mexpt) ,u 2)))
  3482.            ((rat) -1 2))))
  3483.         ,m))
  3484.          (%inverse_jacobi_cn
  3485.           ;; A&S 17.4.52
  3486.           `(($elliptic_f) ((%acos) ,u) ,m)))))
  3487.     (t
  3488.      (recur-apply #'make-elliptic-f e))))
  3489.  
  3490. (defmfun $make_elliptic_f (e)
  3491.   (if (atom e)
  3492.       e
  3493.       (simplify (make-elliptic-f e))))
  3494.  
  3495. (defun make-elliptic-e (e)
  3496.   (cond ((atom e)
  3497.      e)
  3498.     ((eq (caar e) '$elliptic_eu)
  3499.      (destructuring-bind ((fn &rest ops) u m)
  3500.          e
  3501.        (declare (ignore fn ops))
  3502.        `(($elliptic_e) ((%asin) ((%jacobi_sn) ,u ,m)) ,m)))
  3503.     (t
  3504.      (recur-apply #'make-elliptic-e e))))
  3505.  
  3506. (defmfun $make_elliptic_e (e)
  3507.   (if (atom e)
  3508.       e
  3509.       (simplify (make-elliptic-e e))))
  3510.   
  3511.      
  3512.  
  3513. ;; Eu(u,m) = integrate(jacobi_dn(v,m)^2,v,0,u)
  3514. ;;         = integrate(sqrt((1-m*t^2)/(1-t^2)),t,0,jacobi_sn(u,m))
  3515. ;;
  3516. ;; Eu(u,m) = E(am(u),m)
  3517. ;;
  3518. ;; where E(u,m) is elliptic-e above.
  3519. ;;
  3520. ;; Checks.
  3521. ;; Lawden gives the following relationships
  3522. ;;
  3523. ;; E(u+v) = E(u) + E(v) - m*sn(u)*sn(v)*sn(u+v)
  3524. ;; E(u,0) = u, E(u,1) = tanh u
  3525. ;;
  3526. ;; E(i*u,k) = i*sc(u,k')*dn(u,k') - i*E(u,k') + i*u
  3527. ;;
  3528. ;; E(2*i*K') = 2*i*(K'-E')
  3529. ;;
  3530. ;; E(u + 2*i*K') = E(u) + 2*i*(K' - E')
  3531. ;;
  3532. ;; E(u+K) = E(u) + E - k^2*sn(u)*cd(u)
  3533. (defun elliptic-eu (u m)
  3534.   (cond ((realp u)
  3535.      ;; E(u + 2*n*K) = E(u) + 2*n*E
  3536.      (let ((ell-k (elliptic-k m))
  3537.            (ell-e (elliptic-ec m)))
  3538.        (multiple-value-bind (n u-rem)
  3539.            (floor u (* 2 ell-k))
  3540.          ;; 0 <= u-rem < 2*K
  3541.          (+ (* 2 n ell-e)
  3542.         (cond ((>= u-rem ell-k)
  3543.                ;; 0 <= u-rem < K so
  3544.                ;; E(u + K) = E(u) + E - m*sn(u)*cd(u)
  3545.                (let ((u-k (- u ell-k)))
  3546.              (- (+ (elliptic-e (asin (sn u-k m)) m)
  3547.                    ell-e)
  3548.                 (/ (* m (sn u-k m) (cn u-k m))
  3549.                    (dn u-k m)))))
  3550.               (t
  3551.                (elliptic-e (asin (sn u m)) m)))))))
  3552.     ((complexp u)
  3553.      ;; From Lawden:
  3554.      ;;
  3555.      ;; E(u+i*v, m) = E(u,m) -i*E(v,m') + i*v + i*sc(v,m')*dn(v,m')
  3556.      ;;                -i*m*sn(u,m)*sc(v,m')*sn(u+i*v,m)
  3557.      ;;
  3558.      (let ((u-r (realpart u))
  3559.            (u-i (imagpart u))
  3560.            (m1 (- 1 m)))
  3561.        (+ (elliptic-eu u-r m)
  3562.           (* #C(0 1)
  3563.          (- (+ u-i
  3564.                (/ (* (sn u-i m1) (dn u-i m1))
  3565.               (cn u-i m1)))
  3566.             (+ (elliptic-eu u-i m1)
  3567.                (/ (* m (sn u-r m) (sn u-i m1) (sn u m))
  3568.               (cn u-i m1))))))))))
  3569.  
  3570. (defprop $elliptic_eu simp-$elliptic_eu operators)
  3571. (defprop $elliptic_eu
  3572.     ((u m)
  3573.      ((mexpt) ((%jacobi_dn) u m) 2)
  3574.      ;; wrt m
  3575.      )
  3576.   grad)
  3577.  
  3578. (defmfun simp-$elliptic_eu (form y z)
  3579.   (declare (ignore y))
  3580.   (twoargcheck form)
  3581.   (let ((u (simpcheck (cadr form) z))
  3582.     (m (simpcheck (caddr form) z)))
  3583.     (cond ((or (and (floatp u) (floatp m))
  3584.            (and $numer (numberp u) (numberp m)))
  3585.        (let ((u-r ($realpart u))
  3586.          (u-i ($imagpart u))
  3587.          (m (float m 1d0)))
  3588.          (complexify (elliptic-eu (complex u-r u-i) m))))
  3589.       (t
  3590.        (eqtest `(($elliptic_eu) ,u ,m) form)))))
  3591.  
  3592. (defmfun $elliptic_eu (u m)
  3593.   (simplify `(($elliptic_eu) ,(resimplify u) ,(resimplify m))))
  3594.  
  3595. (defun agm (a0 b0 phi)
  3596.   (let (c0)
  3597.     (dotimes (k 16)
  3598.       (psetq a0 (/ (+ a0 b0) 2)
  3599.          b0 (sqrt (* a0 b0)))
  3600.       (setf c0 (/ (- a0 b0) 2))
  3601.       (setf phi (+ phi (lisp:atan (* (/ b0 a0) (lisp:tan phi)))))
  3602.       (format t "~A ~A ~A~%" a0 b0 c0 phi))))
  3603.  
  3604. (defprop %jacobi_am simp-%jacobi_am operators)
  3605.  
  3606. (defmfun $jacobi_am (u m)
  3607.   (simplify `((%jacobi_am) ,(resimplify u) ,(resimplify m))))
  3608.  
  3609. (defmfun simp-%jacobi_am (form yy z)
  3610.   (declare (ignore y))
  3611.   (twoargcheck form)
  3612.   (let ((u (simpcheck (cadr form) z))
  3613.     (m (simpcheck (caddr form) z)))
  3614.     (cond ((or (and (floatp u) (floatp m))
  3615.            (and $numer (numberp u) (numberp m)))
  3616.        ;; Numerically evaluate am
  3617.        (asin (sn (float u 1d0) (float m 1d0))))
  3618.       (t
  3619.        ;; Nothing to do
  3620.        (eqtest (list '(%jacobi_am) u m) form)))))
  3621.   
  3622.